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METHOD AND APPARATUS FOR CREATING A SIMULATED PARTICLE PACK 



BACKGROUND OF THE INVENTION 

Field of the Invention 

[0001] The present invention relates to computerized methods and apparatus for 
creating particle packs and, more specifically, to computerized methods and apparatus for 
creating particle packs with enhanced speed and processing efficiency. 

Description of the Related Art 

[0002] The present invention relates to the creation of particle packs that can be used to 
model or simulate particle containing or particle filled systems. There are many varieties of 
particle containing systems. Specific yet merely illustrative examples are found in the fields of 
asphalts, concretes, ceramics, soils, chemical physics of amorphous materials, slurry mechanics, 
and solid propellants or gas generators. An important class of such systems comprises polymeric 
materials with particles embedded or embodied in the polymer body or binder. 

[0003] It is highly desirable in many circumstances to have a capability to model or 
simulate the properties, performance, etc. of such particle containing systems. In designing such 
particle filled systems, for example, substantial trial and error can be avoided if one is able to use 
a model or simulation to assess the implications of making changes in the components, their 
configuration, the conditions, etc. 

[0004] This is particularly true in microstructural analyses, e.g., wherein the analysis 
focuses on the material at a highly magnified level, such as at the molecular or grain stages. A 
specific example of microstructural analysis would involve the modeling or simulation of a 
particle pack for a propellant, which typically would include fuel and oxidizer particles in a 
polymeric binder. When modeling or simulating a particle containing system, it has become 
customary to construct an analytical replica or representation of the particle filled material, 
which is generally referred to as a particle pack. The term "particle pack" has come to be used in 
the field to refer generally to a collection of particles in a defined volume. 
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[0005] The general object in creating particle packs is to model or simulate the manner 
in which a plurality of particles, usually a large number of such particles, will fill or otherwise 
occupy a particular space. The nature of the particles may vary, in some cases considerably, 
depending on the particular application. These variations may include particle size variations, as 
well as other characteristics. 

[0006] The complexity of models used to simulate real world particle containing 
systems increases dramatically if more than a limited number of variables are considered. The 
numbers of interactions, combinations and permutations brought to bear by anything other than 
relatively simple systems causes the processing requirements to become unwieldy even for the 
most advanced computer systems. This complexity and processing requirement limitation is 
even more pronounced for systems that employ large numbers of particles. 

[0007] These limitations have caused workers in the field to invoke various simplifying 
assumptions to render the problem more workable. One such assumption involves the shape of 
the particle. The complexity and related processing demands can be substantially reduced, for 
example, if one assumes that the particles are spherical. 

[0008] The general approach to simulating particle containing systems in the past has 
involved a more rigorous one in which a particle pack is constructed by placing the individual 
particles into a three dimensional volume, i.e., into a "full field " This approach in theory offers 
promise for accurate simulation of the real particle containing systems. The large number of 
particles required to construct a statistically meaningful or useful particle pack, however, usually 
renders this approach impractical. The processing power required to construct such a particle 
pack is prohibitive for all but extremely small numbers of particles, and only relatively small 
variations in particle size. 

[0009] The simplifying assumptions necessary to reduce processing demands down to a 
manageable level can prevent the resulting particle pack from accurately portraying the actual 
packs that are intended to be simulated. 

[0010] An approach that has been used to make the processing demands more 
manageable and yet the accuracy and reliability of the simulation to be good involves a reduced 
dimension approach. An example of this approach is described in I. Lee Davis and Roger G. 
Carter, "Random Particle Packing by Reduced Dimension Algorithms," J. Appl. Phys., 67(2), 15 
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January 1990 pp. 1022-1029 (hereinafter "Davis and Carter"). The reduced dimension approach 
as described there involves simplifying the pack construction problem by employing the 
principle that an arbitrary straight line drawn through a uniformly random three-dimensional 
particle pack (regardless of whether the particles are spheres) will, on the average, lie inside 
particles the same fraction of its length as the packing fraction. As the length of the line goes to 
infinity, the fraction of its length residing inside the particles of a random pack goes to the 
packing fraction exactly. The line is referred to in Davis and Carter as a "central string." This 
approach has proven quite useful in enabling more complicated particle packs to be constructed. 
It can accommodate a substantially wider range of particle sizes. 

[0011] Notwithstanding the usefulness of the reduced dimension approach to 
constructing particle packs, its direct application has been limited to a certain extent in that its 
accuracy can be compromised by considering only particles intersected by the central string. In 
real three-dimensional packs, and assuming the particles are spherical, each sphere initially 
intersecting the central sphere, as it slides down the central string during placement, comes to 
rest more often on a "nearest-neighbor" sphere than on a lower central string particle. The 
sphere then would roll away from the central string, thereby spreading the spheres over a greater 
distance along the central string. 

[0012] A technique for addressing this limitation, also presented in the Davis and 
Carter paper noted above, involves relaxing the requirement that spheres must intersect the 
central. Under this modified approach, particles are permitted to come to rest within a defined 
range of the central string, after assuming that their path is modified to encounter and avoid 
previously placed spheres. This is referred to generally as "perturbation," or perturbed central 
string packing. 

[0013] This modified reduced dimension approach also has been limited, however, in 
that, when perturbation approaches are used, processing demands once again increase, in some 
cases substantially. This can give rise to a fuller and more complicated particle field that once 
again begins to approach the complexity and processing intensity of a full field. 
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Objects of the Invention 

[0014] Accordingly, an object of the present invention is to provide apparatus and 
methods for creating a particle pack that simulate the actual particle distribution within a space 
with a reasonable degree of accuracy. 

[0015] Another object of the invention is to provide apparatus and methods for creating 
a particle pack that involve reduced processing demands and increased processing efficiency 
relatively to full field simulation. 

[0016] Additional objects and advantages of the invention will be set forth in the 
description that follows, and in part will be apparent from the description, or may be learned by 
practice of the invention. The objects and advantages of the invention may be realized and 
obtained by means of the instrumentalities and combinations pointed out in the appended claims. 

SUMMARY OF THE INVENTION 
[0017] To achieve the foregoing objects, and in accordance with the purposes of the 
invention as embodied and broadly described in this document, a machine-implemented method 
is provided for placing a plurality of particles, each having a characteristic dimension, to create a 
particle pack, the plurality of particles comprising N categories of the particles, the characteristic 
dimension of the particles of a given category being different from the characteristic dimensions 
of the particles of other ones of the categories, the characteristic dimension of the particles 
increasing as the category N increases. The method comprises: a) defining a central string, a 
space disposed about the central string, and N concentric subspaces disposed about the central 
string and within the space, each of the N subspaces corresponding to one of the N particle 
categories; b) selecting a particle from the plurality of particles; c) placing the selected particle in 
the corresponding subspace so that the selected particle becomes a placed particle at a particle 
location unique to that placed particle and is in non-overlapping relation with other placed 
particles, wherein the selected particle placement, according to one aspect of the invention, 
includes defining a catch net representative of buoyancy of a portion of the placed particles and 
positioning the catch net within the space based upon the placement of the portion of the placed 
particles. The selected particle placement according to another aspect of the invention includes 
defining a water level representative of a level of a portion of the placed particles that are smaller 
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than the selected particle and represent a surface of the smaller placed particles, and positioning 
the water level within the space based upon the smaller particle surface. The selected particle is 
placed in non-overlapping relation with respect to the catch net and the water level; and 
d) repeating the particle selection (b) and placement (c) until each of the particles of the plurality 
of particles has become one of the placed particles. 

[0018] In accordance with preferred but optional implementations of the invention, the 
following may apply. The particles comprise spheres and the characteristic dimension of each of 
the particles comprises a radius. The central string comprises a line within the space, preferably 
the z axis of a rectilinear coordinate system imposed within the space. Each of the subspaces has 
a characteristic dimension that corresponds to the characteristic dimension of the particles of the 
corresponding subspace. Each of the subspaces may comprise a cylinder. The subspaces 
comprise concentric cylinders disposed about the central string, each of the cylinders defining 
one of the subspaces and having a radius that corresponds to the radius of the particles in the 
category corresponding to that cylinder. The particle selection may comprise selecting the 
particles randomly or selecting the selected particle from a distribution of the particles. 

[0019] In accordance with another aspect of the invention, an apparatus is provided for 
placing a plurality of particles, each particle having a characteristic dimension, to create a 
particle pack, the plurality of particles comprising N categories of the particles, the characteristic 
dimension of the particles of a given category being different from the characteristic dimensions 
of the particles of other ones of the categories, the characteristic dimension of the particles 
increasing as the category N increases. The apparatus comprises: a) an input device for 
inputting particle selection information; b) a storage device operatively coupled to the input 
device for storing the particle selection information; c) a processor for selecting a particle from 
the plurality of particles, for placing the selected particle at a particle location unique to the 
selected particle in the corresponding subspace so that the selected particle becomes a placed 
particle at a particle location unique to that placed particle and is in non-overlapping relation 
with other placed particles, and for establishing a catch net representative of buoyancy of a 
portion of the placed particles and positioning the catch net within the space based upon the 
placement of the portion of the placed particles, the processor placing the selected particle being 
placed in non-overlapping relation with respect to the catch net. 
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[0020] In accordance with another aspect of the invention, a machine-readable medium 
is provided for use in placing a plurality of particles, each particle having a characteristic 
dimension, to create a particle pack, the plurality of particles comprising N categories of the 
particles, the characteristic dimension of the particles of a given category being different from 
the characteristic dimensions of the particles of other ones of the categories, the characteristic 
dimension of the particles increasing as the category N increase. The machine-readable medium 
comprises: a) machine executable instructions for defining a central string, a space disposed 
about the central string, and N concentric subspaces disposed about the central string and within 
the space, each of the N subspaces corresponding to one of the N particle categories; b) machine 
executable instructions for selecting a particle from the plurality of particles; c) machine 
executable instructions for placing the selected particle at a particle location unique to the 
selected particle in the corresponding subspace so that the selected particle becomes a placed 
particle at a particle location unique to that placed particle and is in non-overlapping relation 
with other placed particles, the selected particle placement instructions including instructions for 
defining a catch net representative of buoyancy of a portion of the placed particles and 
positioning the catch net within the space based upon the placement of the portion of the placed 
particles, the selected particle being placed in non-overlapping relation with respect to the catch 
net; and d) machine executable instructions for repeating the particle selection (b) and placement 
(c) instructions until each of the particles of the plurality of particles has become one of the 
placed particles. 

[0021] In accordance with another aspect of the invention, an apparatus is provided for 
placing a plurality of particles, each particle having a characteristic dimension, to create a 
particle pack, the plurality of particles comprising N categories of the particles, the characteristic 
dimension of the particles of a given category being different from the characteristic dimensions 
of the particles of other ones of the categories, the characteristic dimension of the particles 
increasing as the category N increases. The apparatus comprises: a) an input device for 
inputting particle selection information; b) a storage device operatively coupled to the input 
device for storing the particle selection information; c) a processor for selecting a particle from 
the plurality of particles, for placing the selected particle in the corresponding subspace so that 
the selected particle becomes a placed particle at a particle location unique to that placed particle 
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and is in non-overlapping relation with other placed particles, for establishing a water level 
representative of a level of a portion of the placed particles that are smaller than the selected 
particle and represent a surface of the smaller placed particles, and for positioning the water level 
within the space based upon the smaller particle surface, the selected particle being placed in 
non-overlapping relation with respect to the water level. 

[0022] In accordance with still another aspect of the invention, a machine-readable 
medium is provided for use in placing a plurality of particles, each particle having a 
characteristic dimension, to create a particle pack, the plurality of particles comprising N 
categories of the particles, the characteristic dimension of the particles of a given category being 
different from the characteristic dimensions of the particles of other ones of the categories, the 
characteristic dimension of the particles increasing as the category N increases. The 
machine-readable medium comprises: a) machine executable instructions for defining a central 
string, a space disposed about the central string, and N concentric subspaces disposed about the 
central string and within the space, each of the N subspaces corresponding to one of the N 
particle categories; b) machine executable instructions for selecting a particle from the plurality 
of particles; c) machine executable instructions for placing the selected particle at a particle 
location unique to the selected particle in the corresponding subspace so that the selected particle 
becomes a placed particle at a particle location unique to that placed particle and is in 
non-overlapping relation with other placed particles, the selected particle placement instructions 
including instructions for defining a water level representative of a level of a portion of the 
placed particles that are smaller than the selected particle and represent a surface of the smaller 
placed particles, and positioning the water level within the space based upon the smaller particle 
surface, the selected particle being placed in non-overlapping relation with respect to the water 
level; and d) machine executable instructions for repeating the particle selection (b) and 
placement (c) instructions until each of the particles of the plurality of particles has become one 
of the placed particles. 

BRIEF DESCRIPTION OF THE DRAWINGS 
[0023] The accompanying drawings, which are incorporated in and constitute a part of 
the specification, illustrate presently preferred embodiments and methods of the invention and, 
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together with the general description given above and the detailed description of the preferred 
embodiments and methods given below, serve to explain the principles of the invention. 

[0024] Fig. 1 is a pictorial diagram of a computer system which comprises a presently 
preferred embodiment of the apparatus according to the invention, and upon which the presently 
preferred version of the inventive method may be carried out. 

[0025] Fig. 2 is an illustrative example of a particle pack construction according to the 
preferred implementation, wherein spherical particles are placed along the central string. 

[0026] Fig. 3 shows another illustrative example of a particle pack construction 
according to the preferred implementation, wherein spherical particles are placed along the 
central string, but wherein the particles extend more several particle diameters out from the 
central string. 

[0027] Fig. 4 shows a cylindrical space with cylindrical subspaces therein according to 
the preferred embodiment and the preferred method of the invention. 

[0028] Fig. 5 shows an illustrative particle pack constructed using the preferred 
implementation. 

DETAILED DESCRIPTION OF THE INVENTION 
[0029] Reference will now be made in detail to the presently preferred embodiments 
and methods of the invention as illustrated in the accompanying drawings, in which like 
reference characters designate like or corresponding parts throughout the drawings. It should be 
noted, however, that the invention in its broader aspects is not limited to the specific details, 
representative devices and methods, and illustrative examples shown and described in this 
section in connection with the preferred embodiments and methods. The invention according to 
its various aspects is particularly pointed out and distinctly claimed in the attached claims read in 
view of this specification, and appropriate equivalents. 

[0030] In accordance with one aspect of the invention, an apparatus is provided for 
creating a particle pack. In accordance with another aspect of the invention, a 
machine-implemented method is provided for placing a plurality of particles to create a particle 
pack. 
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[0031] A particle as the term is used herein is used according to its common but broad 
meaning in the field. It includes any type of particle that might be found in a particle rigid body 
of arbitrary shape containing materials. In the preferred implementation, particles comprise rigid 
or nearly rigid objects of arbitrary size and shape. 

[0032] It is useful in analyzing particle packs to be able to segregate the particles into 
categories. A group of substantially identical particles, or a group of particles that for purposes 
of analysis may be deemed identical, and which group is in some way distinguishable from other 
groups of particles is called a "mode." If there is only one type of particle comprising the pack 
and those particles are identical to each other, or may be assumed to be identical, the pack is said 
to be "monomodal." If there are two types of particles, the pack is "bimodal," and if multiple 
types of particles are present, it is "multimodal." 

[0033] A mode represents a principal size or type of particle. In practice, particles even 
within a mode often are not exactly identical, but instead have slight variations. In some 
instances, for example, each mode can be divided into a size distribution about some mean. The 
distribution may comprise a number of discrete sizes about that mean. The sets of various 
particles of a given size or character within a particular mode are called "submodes." In the 
presently preferred implementation, if there is no size distribution attached to a particular mode, 
the submode arrays simply contain the mode data as their first and only entry. Because both 
modes and submodes theoretically merely involve segregating particles into categories based 
upon one or more distinguishing characters, modes and submodes are referred to generally herein 
as "categories." 

[0034] The apparatus according to the invention preferably comprises a computer 
system 100, such as a commercially available personal computer, small business computer, 
engineering work station, mainframe, microprocessor, or virtually any other machine with 
similar processing power, appropriately configured and programmed with computer software or 
"code" as outlined and described herein for performing the tasks identified herein. The computer 
system 100 includes a processor 102, storage in the form of RAM 104, a hard drive 106, a 
compact disk (CD) drive 108, and a drive for diskettes 110, all contained within a cabinet 112. 
System 100 also includes a monitor 116, a keyboard 118, and a pointing device 120 such as a 
mouse, track ball or the like. System 100 also may be part of a network, in which case it may 
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include a network connection 122 and network bus 124. The method according to the present 
invention, in its presently preferred implementations, may be carried out on computer system 
100, using the same hardware and software, as described herein with respect to the preferred 
embodiment of the apparatus according to the invention. 

[0035] "Machine" as the term is used herein refers broadly to a computer, processor, 
microprocessor, a network of such devices, or other apparatus capable of performing processing 
as described herein. Such machines typically and preferably will comprise a computer, such as 
computer system 100. 

OVER VIEW OF THE PREFERRED EMBODIMENT 

[0036] The method and apparatus according to the presently preferred yet merely 
illustrative versions of the invention will be described in this document as involving the 
implementation and execution of a computer program that can be run on a computer system or 
machine as described above. The preferred versions of the apparatus and methods according to 
the invention are referred to collectively in this document as the "preferred implementation." 

[0037] The architecture of the preferred implementation is outlined below. In its 
presently preferred form, it is written in double precision C++ or Fortran 90. This is not, 
however, limiting. The code may, for example, be written in FORTRAN, or other suitable 
programming languages. An example of such a computer system is shown in Fig. 1 . The code 
in its preferred version builds a single particle pack and outputs user-specified pack 
microgeometry and statistics. It is modular and can be disassembled into its component 
subroutines for use elsewhere. 

[0038] In accordance with the apparatus aspect of the invention, an input device is 
provided for inputting particle selection information, and a storage device operatively coupled to 
the input device for storing the particle selection information. The information used by the 
preferred implementation will be described throughout this document. The input device may 
comprise any one or any combination of several input means, including a network download, the 
CD reader, a diskette reader, the keyboard, the tracking device, a tape drive, or any other means 
used for inputting data into a machine capable of implementing these aspects of the invention. 
The storage device also may comprise any suitable storage means, and may include RAM, the 
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hard drive, an optical storage device, storage devices on a machine operably coupled to or 
accessible by computer 1 00, or other suitable storage devices readable by a computer directly or 
with translation. 

[0039] The main program of the preferred implementation comprises seven subroutine 
calls. The names and functions of these subroutines are as follows: 



SUBROUTINE 
NAME 


SUBROUTINE FUNCTION 


INITERROR 


Initializes arrays containing the number of errors, warnings, and 
notations accumulated during program execution. It also initializes 
arrays containing information on how these data are to be 
displayed. 


INNAME 


Looks for and reads an input file named INNAME.INP which 
contains the names of files containing the input data needed by the 
preferred implementation to build the user-specified pack and 
output the requested pack information. If the file INNAME.INP 
does not exist, default file names are used to find instructions on 
pack building. Default values are used if user-specified 
instructions are lacking. If essential instructions are lacking, pack 
building is aborted and error messages are issued. 


SETHISTOGRAM 


Reads the information for the types, quantities, and sizes of 
particles to be used in constructing the pack. 


SETCONTROL 


Reads flags and parameters controlling the construction of the 
pack. Default values are used if files containing these instructions 
cannot be found. 


SETOUTPUT 


Reads flags and parameters controlling the type of pack data to be 
output. Default values are again used if files containing these 
instructions cannot be found. 


PPACK 


Constructs the particle pack according to instructions read by the 
previous subroutines. Each time it is called to build a new pack, it 
initializes all variables and arrays so there will be no interfering 
data left over from a previous call. 


OUTSTAT 


Outputs the user-specified information on pack microgeometry and 
statistics. 



[0040] In summary, the first subroutine (INITERROR) initializes error, warning, and 
notation arrays and parameters. The next four subroutines (INNAME, SETHISTOGRAM, 
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SETCONTROL, SETOUTPUT) and the seventh subroutine (OUTSTAT) control the reading of 
building instructions and the output of data. The sixth routine (PPACK) constructs the particle 
pack. This modularity facilitates incorporation of the code or portions of it into larger calling 
codes, such as codes using numerical methods or finite element analysis. 

[0041] Because in the preferred implementation a substantial amount of information 
may be used to govern particle pack construction, the manner in which the data used by the 
software code is categorized in the context of the preferred implementation will now be 
described. These data categories are placed into common blocks and passed from one code 
section to another as needed. In the presently preferred implementation, all data about the pack 
along with building instructions are placed into one of the following categories: 



DATA BLOCK 


DESCRIPTION 


HISTO 


Governs the types and quantities of particles to be used in 
constructing the pack, i.e., the particle histogram. 


CONTROL 


Contains all parameters and flags controlling how the pack will be 
built and stored. 


CURINFO 


Contains information on the status of the particle currently being 
placed, such as its position and which other particles it is currently 
touching. 


PLACED 


Contains information on all particles placed so far in the pack 
including the position, size, and type of each particle. 


SETOUT 


Contains the flags and parameters governing what pack statistics 
are to be generated and how the pack microgeometry and its 
statistics are to be output. 


ERRORS 


Contains the names errors, warnings, and notations that can occur 
in program input, execution, and output. It is used to monitor the 
number of times each error, warning, or notation occurs and 
records such information in a file and/or to the screen, according to 
the dictates of the user. 



THE SUBROUTINE INITERROR 

[0042] In the presently preferred implementations, there are three categories of 
information, called "advisors," used to monitor the health of the pack, i.e., to monitor whether 
any of the pack construction rules, such as the prohibitions against particle overlap, and general 
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requirements for particle distribution in the pack (e.g., homogeneity at the macro level of the 
random or other specified statistical distributions at the microphone mount level) have been 
violated. They are contained in the program's error arrays and parameters. The three categories 
of advisors are errors, warnings, and notations. 

[0043] Error advisors are of two kinds, fatal and nonfatal. A fatal error will stop 
program execution because there is not enough information available to build the pack. Nonfatal 
errors will not stop the program execution but will use default values in place of erroneous 
user-data. Therefore, a pack will be built but may not be the kind of pack the user intended. 

[0044] Warning advisors are notifications of unusual situations or potential problems. 
If the warning deals with user input, the user should correct the input problem because the 
warning sometimes precipitates the use of defaults that may not have been intended. 

[0045] The third category is notation advisors. There are many data that may be of 
interest during pack building. Notations report on these items. They are not errors or warnings 
but simply notifications of items of interest as pack building proceeds or after it has finished. 

[0046] All advisors in this implementation have arrays that keep track of how many 
times each has been called: NumErrors(z), NumWarnings(z), and NumNotes(z), for the z-th error, 
warning, or notation, respectively. The totals for each category of advisor are also accumulated 
in the parameters NTotErrors, NTot Warnings, NtotNotes. 

[0047] The display of advisors in this preferred implementation is controlled by arrays 
containing integers limiting the number of times an advisor is to be displayed to the computer 
screen and written to a disk file created for that purpose. The integer arrays LimScErrors(z) and 
LimDiErrors(z) dictate how many times the z-th error is to be written to the screen and the disk 
file, respectively. Similarly, LimScWarnings(z) and LimDiWarnings(z) are provided for 
warnings and LimScNotes(z) and LimDiNotes(z) for notations. If any of these array elements is 
less than or equal to zero, the corresponding advisor is not written. If the user wants all advisors 
written every time they are encountered, the limit arrays need simply be set to a very large 
integer (in this implementation not to exceed 2147483646, which is one less than the maximum 
integer size). The total number of times any advisor has been encountered is monitored in 
NTotErrors, NTotWarnings, and NTotNotes whether or not it is displayed. 
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[0048] All data on all advisors according to this preferred implementation are kept in 
and monitored by the subroutine ErrWarnNote(zerr ? iwam, inote, iunitnumber), where ierr is the 
integer specifying the error, and similarly for iwarn and inote. The integer iunitnumber is the 
disk file unit number to which advisors are written as requested. It is assigned the number 13 in 
this implementation. 

THE SUBROUTINE INNAME 

[0049] In the presently preferred implementation, the data required by the code are read 
from several data files. The INNAME subroutine reads the names of the data files in which the 
data resides from a file called INNAME.INP. The files whose names are specified in 
INNAME.INP are as follows: 



FILE 


DESCRIPTION 


DEFAULT NAME 


Histogram 


Contains data on sizes, quantities, and types of particles 
along with material information and the sizes of the 
cylinders into which they are to be dropped 


HIST.INP 


Drop File 


Contains data on where and how particles are to be 
dropped onto the pack. 


DROPDATA.INP 


Preplace File 


Containing the positions and types of particles to be 
preplaced in the pack, if any. 


PREPLACE.INP 


Input Control 
File 


Contains the input flags and parameters that control how 
the pack is to be built. 


INFLAG.INP 


Output Control 
File 


Contains flags and parameters controlling what output 
data is to be generated and how it is to be displayed. 


SETOUT.INP 



[0050] In accordance with the preferred implementation, if the file INNAME.INP does 
not exist, then default names are used for each of the five files containing input data. The default 
names for this implementation are identified in the third column of the table above. 

[0051] In the preferred implementation, the histogram file is mandatory but the other 
files are not. Packs can be built if all other files are missing by using defaults. Of course, the 
various options and data production capabilities of the presently preferred implementation may 
not be exercisable without these other files. 
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THE SUBROUTINE SETHISTOGRAM 

[0052] The SETHISTOGRAM subroutine according to this preferred implementation 
reads the particle histogram from the histogram file and sets necessary parameters in the program 
associated with that file. The particle histogram comprises a set of particle modes and 
submodes, in this implementation in the form of arrays. If there is no size distribution attached 
to a particular mode, the submode arrays simply contain the mode data as their first and only 
entry. 

[0053] First, subroutine SETHISTOGRAM checks to see if the histogram file exists. If 
it does not, the program is terminated with a fatal error message. If the file exists, the data is 
read as follows. 

[0054] The number on the first line of the histogram file is the "stop build" criterion in 
the form of an integer called STOPBUILD. If STOPBUILD is positive, then that is the total 
number of particles to be placed in the pack. If STOPBUILD is less than or equal to zero, then 
pack building is to be stopped when the center of the highest placed particle reaches a certain 
altitude or when a prespecified maximum number of particles is placed, whichever comes first. 
If STOPBUILD is less than or equal to zero, the second line is the real number 
STOP ALTITUDE specifying the altitude at which pack building is to cease. 

[0055] The particle histogram is processed next in this implementation. Each set of 
lines contains the data for one mode. The modes do not have to be in any particular order. The 
first line of the set contains the mode information. If there is a size distribution associated with 
this mode, a set of lines follows that specifies the submodes. In this implementation there are six 
numbers on each mode line. 

[0056] 1. The first number is the mean particle size for the mode. It is read 
in as a particle diameter but is internally converted to particle radius. 

[0057] 2. The second number is the mass quantity for the mode. Because the 
total mass quantities are normalized later, relative mass quantities are adequate. 

[0058] 3. The third number is the mass density for the mode. The units of 
particle size and density preferably are consistent. In this implementation this is a requirement. 
For example, one should not input a mode diameter in microns and its density in grams/cm 3 . 
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[0059] 4. The fourth number is an integer specifying a material number for 
the mode. Material numbers are a way of specifying some property unique to a mode or set of 
modes. For example, a material number may specify a particular material composition or it may 
specify a set of modes that have some other feature in common. 

[0060] 5. The fifth number is the radius of the cylinder into which the 
particles of this mode will be dropped. The centers of the particles are constrained to stay within 
this cylinder wall. If there is a distribution associated with this mode, the cylinder radius of the 
mode is ignored and the cylinder radii for each submode are used. 

[0061] 6. The last number is a logical parameter RDIST specifying whether 
there is a distribution of submodes associated with this mode. If the parameter is FALSE, there 
is no distribution and the first submode array elements are filled with information for that mode. 
If the parameter is TRUE, the program next reads the submode information in the following 
form. 

[0062] The submode data for this implementation comprises a set of lines, each one of 
which includes the following three pieces of information: 

[0063] 1. The submode size. As with the mode sizes, it is read in as a 
diameter but internally converted to a radius. 

[0064] 2. The fractional mass of the submode. This number is the fraction of 
the mode's mass contained within this submode. The sum of a mode's fractional masses should 
add to unity. The program normalizes them if they do not. 

[0065] 3. The submode cylinder radius. 
[0066] In many parts of the code as presently preferred and implemented, it is 
convenient to have the submodes ordered in increasing particle size. A global size array is 
created that contains all submodes in increasing particle size. The global submode arrays 
comprise an array of submode radii (QSIZE), submode mass (QMASS), submode mass density 
(DENSITY), submode material number (MODEM AT), and submode cylinder radii 
(SUBCYLRAD). 

[0067] The end of a submode is signaled by a set of zeroes everywhere that a nonzero 
number was expected, and the end of a mode is likewise signaled when a set of zeroes is read. 



17 



[0068] After reading all mode and submode information, the SETHISTOGRAM 
subroutine reads the information in the drop file if it exists. If it does not exist, the program 
drops the particles randomly anywhere within the cylinder corresponding to that mode or 
submode. If it does exist, then it contains the following data. 

[0069] The first line has two integers. The first integer is a drop option, either 1, 2, 
or 3. The second integer is the number of particles that are to be dropped according to that 
option. 

[0070] Option 1 means that first group of particles are to be dropped randomly 
anywhere within their cylinder. 

[0071] Option 2 means that the particle group is to be dropped within a circle of 
relative radius a r where the drop radius is some fraction of the cylinder radius (0 < a r < 1) of the 
particle that will be selected. The center of this circular region is at relative radial position p r 
where p r is some fraction of the particle's cylinder radius when it is chosen (0 <p r < 1). Because 
dropping particles outside the cylinder wall is forbidden, it should hold that p r + a r < 1 . These 
data are stored in the pack position arrays. 

[0072] Option 3 means that each particle in that group is to be dropped from a specified 
position and is of a specified mode. The x position of the drop site is stored in the x position 
array, the y position in the y position array, and the mode number in the mode number pack 
array. If an (x, y) position lies outside the cylinder corresponding to that mode, a warning 
message is issued and the particle is randomly dropped anywhere within the cylinder. 

[0073] The end of this input file is signaled by a value of zero for the drop option. 

[0074] Another function of the SETHISTOGRAM subroutine in this implementation is 
to create an index array that orders the global cylinder radii in ascending order with an integer 
array called CylSizeOrder so that cylrad(CylSizeOrder(l)) is the smallest cylinder radius, 
cylrad(CylSizeOrder(2)) is the second smallest, and so forth. 

THE SUBROUTINE SETCONTROL 

[0075] The currently preferred implementation includes a variety of options for pack 
construction. Some of those options are read in a file that can have any name as specified in 
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Subroutine INNAME but whose default value is INFLAG.INP. This subroutine SETCONTROL 
reads those control parameters. 

[0076] If a data file containing input control parameters exists, then the control 
parameters are read from it in this implementation. If it does not exist, a default is set for each 
parameter in this routine. The parameters used in this implementation are described in detail 
below but in summary include the following: 

[0077] ROUGHNET = A flag that specifies whether the k-th mode catch net (described 
below) is ratcheted up after each particle placement so that no two particles are at exactly the 
same altitude. It is mostly used at the beginning of pack building so the spheres at the bottom of 
the pack do not lie exactly on the bottom of the cylinder. The default is TRUE. 

[0078] POOL = A flag that specifies whether the number of particles in the pack are 
selected by a random number generator, weighted by their frequency of occurrence in an infinite 
pack, or whether they are selected from a pool whose mode proportions are, to the nearest 
integer, exactly the proportions of an infinite pack. If the POOL flag is TRUE, then random 
chance variations in modes that contribute few particles to the pack are avoided. The default is 
TRUE. 

[0079] FINDLOW = A flag that orders find the lowest exposed north pole of placed 
particles within its drop cylinder and to drop immediately in its neighborhood. It is overridden 
by any instructions that may appear in the drop file. The default is TRUE. 

[0080] NDROP = A positive integer that specifies the number of times a trial drop of 
each particle is to be made. The lowest final resting place is selected from among the trial drops. 
Higher values for NDROP means the pack will be more dense. The default is one. 

[0081] ISEED = An integer specifying the random number seed. A non-positive 
integer will trigger to seed the random number generator from the system clock. 

[0082] TUNNEL = A flag that gives permission to check for cavities underneath larger 
spheres when a smaller sphere is dropped. The small sphere "tunnels" through the large one and, 
if there are no overlaps in the southern or underneath side, is dropped again from that underneath 
position whereupon it searches for a final resting place in that cavity. The tunneling procedure 
seeks to mimic tamping or shaking a pack to fill in such cavities, thereby making the pack more 
dense. 
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THE SUBROUTINE SETOUTPUT 

[0083] A particle pack is typically created so that the system or method user can obtain 
specific information about the actual particle pack the simulation is intended to simulate. The 
information required depends on the application. The subroutine SETOUTPUT in this 
implementation reads the flags and parameters that instruct other parts of the preferred 
implementation to create specific pack information. In this implementation, for example, if the 
instruction file for setting output flags and parameters does not exist, the code outputs a limited 
amount of information into a file called PACKDATA.OUT. Otherwise, this basic pack 
information is output to a file name specified in Subroutine INNAME. With no specifications, 
the following data is sent to the output file: 

• an echo of the particle histogram; 

• the packing fraction; 

• the maximum and minimum positions of particles for each submode as 
well as overall maxima and minima; and 

• the number of particles placed for each submode as well as total particles 
placed. 

[0084] If the pack output file exists when using the preferred implementation, it 
preferably contains the following information. The first line contains a flag specifying whether 
the particle pack is to be written to a disk file. The default is FALSE. If the flag is TRUE, then 
the next set of lines specifies how the pack file is to be written. On the first subsequent line, the 
format of the pack (ASCII or binary) is written. On the second subsequent line, a flag is read 
that determines whether the material number and radius are written in the pack file in addition to 
the three coordinates and the mode number. On the third subsequent line, the output geometry 
(Cartesian or cylindrical coordinates) is written. The name of the file to which the pack is to be 
written is on the fourth line. 

[0085] The next line contains a flag specifying whether radial distribution functions 
gij{r) are to be generated. If it is TRUE, then the functions that are to be generated are specified 
in a list of integer pairs i and j on the subsequent lines, where the / and j represent modes or 
submodes. Each line also contains a Ar which is the increment in terms of the y-th particle radius 
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in which the number of particle centers are to be counted (default = 0.1 radius) and r max which is 
the extent to which the radial distribution function extends in terms of y'-th particle radius. Hence 
each radial distribution function requested is on a line with the format 

i j Ar r max . 
The list ends when zeroes are read. Each radial distribution function is written to a file with the 
name #[*]_[/]. dat where each set of square brackets indicates an actual specified integer to be 
used in the file name, e.g., g3_14.dat for the radial distribution function of Mode 14 spheres 
about Mode 3 spheres. The default for the radial distribution flag in this implementation is 
FALSE. 

[0086] The next line in this implementation contains a flag specifying whether 
coordination numbers are to be generated. If it is TRUE, then the coordination numbers that are 
to be generated are read as a list of integer pairs i and j on subsequent lines. The list ends when 
zeroes are read. The coordination numbers preferably are all written to one disk file 
COORNUM.DAT. The default for the coordination number flag in this implementation is 
FALSE. 

[0087] Radial distribution functions and coordination numbers and their calculation 
procedures are defined in the discussion of subroutine OUTSTAT, below. 

[0088] In accordance with the apparatus aspect of the invention, the apparatus 
comprises means for defining a space, a plurality of subspaces within the space, and a central 
string within the space and within the subspace. In accordance with the method aspect of the 
invention, the method comprises defining a space, defining a plurality of subspaces within the 
space, and defining a central string within the space and within the subspaces. 

[0089] The subspaces preferably comprise N concentric cylinders disposed about the 
central string and within the space. Each of the N cylinders corresponds to one of the N particle 
modes and has a cylinder radius corresponding to the particle radius of the particles of the 
corresponding particle mode. The means as implemented in the preferred embodiment 
comprises the processor. 

[0090] The space and subspaces are constructions that define the regions or volumes 
into which the respective particles are to be placed. The space may comprise any contained 
geometric region, but preferably comprises a cylindrical construction. The subspaces also may 



21 



comprise any contained geometric region, but preferably have a shape and configuration that 
corresponds to the overall space. In the preferred versions of the apparatus and method, the 
space comprises a cylindrical volume, and the subspaces comprise a plurality of concentric 
cylinders, disposed about a longitudinal or z axis, wherein sequential ones of the cylinders have a 
radius that is larger than the prior cylinders. An illustrative example of such a concentric 
cylinder configuration is shown in Fig. 2. The number of cylinders in the preferred 
implementation is equal to the number of particle categories, e.g., submodes, that will be present 
in the particle pack. The radius of each cylinder is selected to correspond to the radii of the 
particles. The radius of the first cylinder corresponds to or is a function of the particle radius of 
the smallest mode, the radius of the second cylinder is equal to the particle radius of the second 
smallest mode, and so on. The cylinder radius for a particular particle category preferably 
corresponds to or is a function of the radius or a representative radius of the corresponding 
particle for that mode or submode. The cylinder radius may be any positive number multiple of 
the particle radius. Preferably, but optionally, the cylinder radius is about 3 to 10 times the 
particle radius, and more preferably, about 5 to 10 times the particle radius. 

THE SUBROUTINE PPACK 

[0091] Further in accordance with the method aspect of the invention, the method 
includes selecting a particle from the plurality of particles, and placing the selected particle in the 
corresponding subspace so that the selected particle becomes a placed particle at a particle 
location unique to that placed particle and is in non-overlapping relation with other placed 
particles. 

[0092] In accordance with the preferred implementation, particles are selected 
randomly from the distribution of particles using a Monte Carlo method, although this is not 
necessarily limiting. In terms of the apparatus, this comprises a random number generator in 
computer system 100. 

[0093] Regarding terminology, a "selected particle," also referred to herein as a 
"current particle," is a particle that has been selected for placement, but has not yet been placed 
(become a "placed particle") in a particle location in the particle pack. A "placed particle" is a 
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particle that has been placed or fixed by assigning it a unique particle location in the particle 
pack. 

[0094] In the preferred implementation, the software code subroutine selects particles 
and constructs the particle pack according to instructions read by the previous subroutines. Each 
time it is called to build a new pack, it initializes all variables and arrays so there will be no 
interfering data left over from a previous call. Code lines of PPACK and a brief explanation of 
each will now be provided. 

CALL PREPAREPACK 

[0095] Subroutine PREPAREPACK performs three functions in this implementation. 
First, it preplaces particles, if any, that have a position assigned by the user, such as would be 
desirable if semirandom packs were to be built. Secondly, it arranges histogram data such that 
particles can be chosen by the random number generator with probabilities appropriate for their 
abundance in the reduced dimension pack. Third, it initializes and sets parameters that will be 
used by the rest of the subroutines in PPACK. 

DO icarus = 1, stopbuild 

[0096] This DO loop continues to cycle until either the user-specified number of 
particles has been placed in the pack or until a pack height has been reached. 

CALL CHOOSE 

[0097] Subroutine CHOOSE chooses the type, category, mode or submode of the next 
particle to be dropped. It drops the "current particle" according to user-specified instructions. 

DO idrop = 1, ndrop 

[0098] The user may specify more than one drop at different drop sites for each particle 
so that the lowest final site can be selected. Multiple drops are equivalent to tamping the pack so 
the particles settle into lower sites resulting in a denser pack. 
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CALLDROPSITE 

[0099] Subroutine DROPSITE selects an (X, Y) position from which the particle will be 
dropped. The site can be chosen randomly, or the user may ask the code to use a criterion to find 
potential sites where the particle is likely to settle to an especially low final resting place. 

CALL PLACE 

[00100] Subroutine PLACE finds the final resting place or "placed particle location" for 
the current particle, which is the site at which the current particle will ultimately be placed. The 
subroutine PLACE in this implementation contains all the instructions on how a particle drops, 
rolls, and encounters other objects in the pack as it moves toward its final resting place. It is 
where the machine spends the great majority of its processing time during pack construction. 

END DO 

CALL UPDATE 

[00101] Subroutine UPDATE updates all the particle pack arrays and parameters with 
the position and additional data associated with the newly placed particle. 

END DO 
RETURN 
END 

[00102] In accordance with the preferred implementation, a particle pack is constructed 
by "dropping" the selected particles one at a time into the set of concentric cylinders. The 
dropping presumes a vertical orientation, which is not necessarily required. But regardless of the 
reference frame, the current particle is moved in the space until it encounters another object, as 
described in principle herein. Each particle mode is dropped somewhere within the cylinder 
radius of its corresponding cylinder. The radius of the cylinders in this preferred implementation 
is specified by the user, but preferably the small particles are dropped into a small-radius 
cylinder and the larger particles into a larger-radius cylinder or subspace. The large particles can 
therefore overlap the inner small-radius cylinders but the small particles do not have to fill the 
volume in the larger cylinders. The reduced dimension approach as described in Davis and 
Carter is used in this preferred implementation. 
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[00103] The reduced dimensional approach can make pack building tractable from a 
processing demand perspective when the particles comprising the pack have broad size ranges. 
Particle centers in this implementation cannot pass outside the radius of the cylinder wall to 
which they belong but can pass through any smaller cylinder walls. Because only the innermost 
cylinder can have all modes or submodes intersecting it, only there does the reduced dimension 
pack literally approximate the corresponding real three-dimensional pack. The remainder of the 
pack can however be simulated using this approach, for example, by assuming that the remainder 
of the space is statistically identically configured. 

[00104] In the presently preferred implementation, the selection of particles to be placed 
and the placement of the particles in the particle pack are undertaken by the processor within the 
machine under the execution of the PPACK subroutine. In this implementation, there are four 
subroutines called by PPACK (PREPAREPACK, CHOOSE, PLACE, and UPDATE). Each of 
these subroutines and the functions they perform will now be described. 

THE SUBROUTINE PREPAREPACK 

[00105] As noted above, the PREPAREPACK subroutine has three functions: 
preplacing particles, setting cumulative arrays so the random number generator can readily 
choose the mode of the next particle to be placed, and initializing various parameters in 
preparation for building the pack. Each of these functions will now be described in turn. 

Preplacing Particles 

[00106] The first function of PREPAREPACK is to preplace any particles for which the 
user has specified a position. This option permits semirandom or fully ordered packs to be 
constructed. A filename is read in Subroutine INNAME that contains the global (X 9 7, Z) 
position and the mode or submode number of each particle to be preplaced. If the filename is not 
read, the default filename for preplaced particles is PREPLACE.INP. If a filename of preplaced 
particles does not exist, the program assumes no particles are to be preplaced and this option is 
not exercised. The preplace file contains five numbers for each particle to be placed. In order, 
they are the X position, the Y position, the Z position, and the mode number, and the submode 
number. These are read into the surface arrays. 



25 



Setting Arrays for Choosing Particles 

[00107] The second function of PREPAREPACK according to the preferred 
implementation is to make an array so that the random number generator can randomly choose, 
with the appropriate probability, the submode of the next particle to be placed in or on the pack. 
The relative fraction of particle numbers for each submode therefore should be calculated. The 
number of particles in each submode to be placed will depend on the submode mass, its mass 
density, its particle radius, and its cylinder radius, e.g., in the following manner. 

[00108] In this implementation there are two options in setting up the weighting factors 
that govern how the random number generator will choose the next particle to be placed. The 
first option is to use the same weighting factors each time a particle is chosen for placement. 
This is referred to herein as the "unchanged-weighting option." The second option, referred to 
herein as a "pool option," creates a pool that has a set of bins with the appropriate number of 
particles placed into it. Each time a particle is drawn from one of the bins, the weighting factor 
of drawing from each bin is altered slightly. Creating a pool of particles from which to draw 
forces the right proportion (to within nearest integers) of particles from each submode to be 
placed. In the limit of an infinite number of particles in the pack, these two methods are 
identical. The pool option more closely resembles a very large pack. The unchanged-weighting 
option will be described first. 



Unchanged- Weighting Pool Option 

[00109] The total volume of k-th submode particles to be placed is 



(1) 



where m k is the submode mass, pk is the submode mass density, and a* is the particle radius for 
that submode. 

[00110] The relative number of submode k particles in this reduced dimension pack is 
also proportional to the area of the confining cylinder for that submode: nR k 2 where R k is the 
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radius of the cylinder. Therefore, the fraction f k of &-th submode spheres to be placed in a 
reduced dimension pack is proportional to 



m k R 2 k 



(2) 



[001 11] To obtain the actual fraction the above expression can be normalized: 



F=£<t> k (3) 



where N is the total number of particle submodes. Then 

f k =<pJF. (4) 

These are the unchanged-weighting factors that determine which particle submode shall be 
chosen for the next placement. 

[00112] These weighting factors are placed into a cumulative array so the random 
number generator can immediately choose the submode. In other words, the cumulative array is 
the sum of all weighting factors prior to it: 

C, = /, 

Q=tf i9 \<i<N (5) 
k=i 

C N = 1. 

[00113] When a uniform random number r is generated on [0, 1], a bisection search is 
made of the array {C/, C2, . . . C^to find which pair of numbers the random number lies 
between. The A:-th submode is chosen if C*./ < r < C*. Co is defined to be zero. 
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Pool Option 

[00114] The pool option is similar. As discussed before, the user can specify the 
number of particles to be placed in the pack or can specify a pack height in which case the 
program continues to add particles until the height is achieved (or until the array sizes are 
exhausted). The pool option in this implementation can only be used if the user has specified the 
total number of particles to be placed in the pack, say n. The number of submode k particles n k 
to be placed is 

n k = mrr[nf k ] (6) 

where NINT indicates the integer nearest to its argument. Because of nearest integer roundoff, it 
is possible that the n k do not sum exactly to n. Then the quantities should be adjusted to 
compensate where the roundoff errors were the largest. 

[00115] It is possible, especially for packs with very broad particle size distributions, 
that some nfk < Vi so that submode would not be represented at all in the pack. The preferred 
implementation issues a warning when that happens, or even when a particular is very small, 
so that the user is aware that the pack is too small to be able to sample enough different 
configurations to mimic a true three-dimensional particle pack. Changing cylinder sizes on some 
submodes can solve this situation, should it arise. 

Initializing Various Parameters 

[00116] The packing routine can be called many times, as when the user wishes to 
search for an optimum pack subject to some set of constraints. Some of the parameters and 
arrays from previous calls must not be carried into successive calls or they can taint the results. 
This portion of Subroutine PREPAREPACK initializes all necessary parameters, flags, and 
arrays for use in a new pack simulation. 

THE SUBROUTINE CHOOSE 

[00117] The CHOOSE subroutine uses the random number generator to choose the 
submode of the next particle to be added to the pack. The random number generator in this 
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implementation generates a uniform distribution on [0, 1], The cumulative array from the 
previous routine PREPAREPACK is a set of numbers such that 

Co = 0 < C } < C 2 < ... C N -i <C N =l. (7) 

[00118] Because this cumulative array can be large, and because there is no restriction 
on the relative distances between successive array elements, in this implementation, a search by 
bisection is a preferred method to find two numbers between which a given random number lies. 

[00119] When a random number r is chosen in [0, 1] we want to find the two array 
elements that bracket it such that 

C k .j<r<C k A<k<N. (8) 

(In the unlikely event that r = 0, it will be assumed to belong to the first bin, />., where k = 1.) 
In preferred bisection methods, a check is made to determine on which side of Cn/2 the random 
number lies. Then attention is concentrated on that side and it is bisected again, continually 
narrowing the range of {Ck} in which r lies until Eq. 8 is satisfied. 

[00120] As noted above, the particle just chosen is referred to herein as the "current 
particle" or the "selected particle" in the discussion below. 

THE SUBROUTINE DROPSITE 

[00121] The DROPSITE subroutine chooses a drop site from which the current particle 
can be dropped. Actually, if the user so specifies, it chooses several drop sites and then keeps 
the site at which the particle came to the lowest possible final resting place, referred to herein as 
a "placed particle location." In a sense, this is the same as tamping the pack so that it settles into 
a denser configuration. The user specifies how many trial drops a particle is to undergo. The 
default in this implementation is one. 

[00122] There are several options available to the user in this implementation when 
deciding where a new selected particle is to be dropped onto the pack. 

[001 23] 1 . Drop it randomly anywhere within its cylinder; 
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[00124] 2. Drop it randomly anywhere within a user-specified circle that is 

contained within the cylinder; 

[00125] 3. Drop it from a user-specified (X, Y) position; and 

[00126] 4. Select the lowest exposed north pole among the placed pack 

spheres and drop it randomly within a short distance of that pack sphere's (X, Y) position. If 

dropping multiple times, it is preferable to select the several lowest exposed north poles and drop 

within a short distance of each of them. 

Dropping Randomly Within the Cylinder 

[00127] When dropping a sphere anywhere within a given cylinder, two random 
numbers are generated, u 9 and u p . The azimuthal angle (in radians) of the drop position <p c 
relative to the central string is chosen by 



[00128] The radial position is not chosen uniformly because there is more area for large 
radial positions in a circle than for small radial position. Because the area grows as the square of 
the radial size, all areas will be equally accessible if we take the square root of the random 
number generator and multiply it by the cylinder wall radius for the current sphere W c \ 



(p c = 2nu 9 . 



(9) 




(10) 



[00129] The (X, Y) drop site is found as follows: 



X c = p c cos (p c 



(11a) 



Y c = p c sin <p c . 



(lib) 
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Dropping Randomly Within a User-Specified Circle 

[00130] Suppose the user specifies that the particles are to be dropped randomly within a 
circle of radius R centered at (Xo, Jo). In the preferred implementation it is required that the 
circle be contained within the cylinder of the current sphere, i.e., 



J X\ + Y\ +R<W C . (12) 

Otherwise, the circle is an illegal specification because a particle cannot be dropped outside its 
cylinder wall. In that case, the preferred implementation will default to dropping randomly 
anywhere within the cylinder wall at W c . A value of zero for R is permitted, meaning the 
particles are dropped from the position (Xo, Y 0 ). If this condition is satisfied, then the drop site 
can be selected as follows. Two random numbers are generated, u 9 and u p . Then 

(p c = 2%u 9 (13) 

Pc = R^ P (14) 
X c =X 0 +p c cos <p c (15a) 



and 



Y C = Y 0 + p c sin <p c . (15b) 
Dropping From a User-Specified AT Position 

[00131] The condition of dropping a selected particle from a user-specified X-Y position 
is a subset of the previous discussion for the special case of R = 0. The positions are stored in the 
pack surface arrays for that particle until it is placed, at which time the drop XY position is 
replaced by the XY position of the final resting location, or the placed particle location, for that 
particle. 
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Dropping Above Lowest Sites 

[00132] In the preferred implementation, a running list is kept of all placed pack spheres 
whose north poles have not yet been covered by other pack spheres. These north poles are 
ordered according to their altitude (north pole altitude, not sphere center altitude). It is generally 
desirable for the current or selected sphere to find a final resting place with the lowest altitude. 
The lowest exposed north pole is a more likely candidate for a low final placed particle location. 
Therefore, if the current sphere is to be dropped m times (m > 1), then we choose the lowest m 
exposed north poles and drop randomly within a circle of radius MIN(a c , ai) of that north pole's 
XY position, as long as it is within the current sphere's corresponding cylinder wall. The 
quantity a c is the radius of the current sphere and is the radius of the placed sphere above 
whose exposed north pole the current sphere is being dropped. Dropping randomly within this 
circle is done in the same manner as is described above. North poles outside W c are ineligible for 
dropping in this implementation. If there are less than m exposed north poles, the other sites are 
chosen randomly within the cylinder. 

THE SUBROUTINE PLACE 

[00133] The subroutine PLACE in this implementation takes the current particle from its 
initial drop site to the final resting place or placed particle location associated with that drop site. 
The dropping and rolling preferably is done without any dynamics. The particles in this 
implementation do not bounce, nor do they gain speed or shoot off the edge of another particle 
with a parabolic trajectory. It is as if the particles are being dropped in a highly viscous fluid so 
that all inertia is absent and the particle creeps to its final resting place. 

[00134] There are four subroutines in the presently preferred implementation of the 
PLACE subroutine. They are referred to herein as FOBJ1, FOBJ2, FOBJ3 and STABLE. 
FOBJ1 finds the first object the descending current sphere encounters. FOJB2 finds the second. 
FOJB3 finds the third. STABLE determines which contacted objects the current sphere will stay 
in contact with and which it will roll away from. The term "object" is used instead of "sphere" 
because one of the objects of contact may be the cylinder base or wall. 

[00135] A number of terms will now be defined for purposes of this discussion before 
describing the operation of the PLACE subroutine. The "roll corridor" is the path of descent of a 
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sphere as it descends along its contacts with one or more spheres and/or the cylinder wall. A 
good share of the effort of subroutine PLACE is to find the objects that overlap the roll corridor. 
An "object" according to this implementation can be one of four things that the current sphere 
can encounter as it descends: (1) other pack spheres (placed particles), (2) the cylinder base or 
wall, (3) a "water level," or (4) a "catch net." The latter two are described more fully below. 
"Contact stability" refers to whether the current sphere is in compressive contact with an object 
or in tensile contact. If it is in compressive contact with another sphere or the cylinder wall, they 
are pushing against each other due to the weight of the current sphere when gravity is acting 
along the -z axis. The current sphere at this position stays in contact with an object it 
compressively touches. Usually if it touches three or more objects compressively, it is stable and 
is placed at that position. If there is a tensile contact, the object and current sphere are pulling 
away from each other. Therefore, the current sphere will roll away from the object unless it 
already has three or more compressive contacts. The first pack sphere touched by the current 
sphere is denoted as Sphere 1, the second as Sphere 2, and so forth. If a set or plurality of 
spheres is touched simultaneously, the numerical order can be arbitrarily selected. 

[00136] Subroutine FOBJ1 searches for the first object of contact. If the current sphere 
simultaneously contacts more than one object, it checks them for stability through subroutine 
STABLE. Subroutine FOBJ2 searches for the second object of contact. If the current sphere 
simultaneously contacts two or more objects, it checks them for stability. Likewise, subroutine 
FOJB3 searches for the third object of contact. Again, if the current sphere simultaneously 
contacts three or more objects, it checks them for stability. The current sphere is placed if it is 
stably touching three or more objects or has hit the catch net or a water level. 

[00137] We will now give an overview of each subroutine and how they interact, 
followed by a detailed discussion of each one. 

Overview of Subroutine STABLE 

[00138] Subroutine STABLE determines whether the current sphere is in compressive or 
tensile contact with all objects it is currently touching. It marks those objects in tensile contact 
so that the current sphere will roll away from them as it continues to descend. 
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Overview of Subroutine FOBJ1 

[00139] The first subroutine for determining object contact is called FOBJ1 for "Find 
Object 1," object 1 being the first object the descending current sphere encounters. There are 
four possible outcomes of FOB J 1 : 

[00140] • 1 . The current sphere encounters the catch net or a water level 
before encountering any other sphere. If it does, then it is placed at the position where such 
contact occurs. 

[00141] • 2. The current sphere encounters a placed sphere of the pack. 
If it does, the index number of that sphere (hereafter denoted as Sphere 1) is recorded and FOBJ2 
is called. 

[00142] • 3. The current sphere simultaneously hits two placed spheres 
of the pack. If it compressively touches both of them, it will roll along the roll corridor defined 
by them and FOBJ3 is called to determine whether a third object is hit before it loses contact 
with one or both of the touched placed spheres. 

[00143] • 4. The current sphere simultaneously hits three or more placed 
spheres of the pack. STABLE checks to see if any three of them provide a stable placement. If 
they do, the current sphere is placed there. If they do not, STABLE marks the objects from 
which the current sphere will roll away. 

Overview of Subroutine FOBJ2 

[00144] FOBJ2 searches for the second object to be encountered by the current sphere as 
it descends along the path of steepest descent on Sphere 1. There are four possible outcomes of 
FOJB2: 

[00145] • 1. The current sphere hits the catch net or a water level before 
encountering any other object. If it does, it is placed there. 

[00146] • 2. The current sphere loses contact with Sphere 1 by reaching 
its equator before encountering any other object. If it does, it is free falling again and FOBJ1 is 
called, and re-initialized. 

[00147] • 3. The current sphere hits a second object, either another 
placed sphere or its cylinder wall. STABLE is then called to see if it is in compressive contact 
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with both objects. If it is, FOBJ3 is called. If it loses contact with one object, FOBJ2 is recalled 
and re-initialized. If it loses contact with both objects, it is in free fall again and FOBJ1 is 
recalled and re-initialized. 

[00148] • 4. The current sphere hits two or more additional objects 
simultaneously. STABLE is called to check for compressive contacts. If there are three or more 
compressive contacts, the particle is placed where it is. If there are two, FOBJ3 is recalled; if 
there is one, FOBJ2 is recalled; and if there are none, the current sphere is in free fall and FOBJ1 
is recalled. 

Overview of Subroutine FOB J3 

[00149] FOBJ3 searches for the third object the current sphere encounters as it descends 
along the roll corridor formed by the first two objects. There are four possible outcomes of 
FOBJ3: 

[00150] • 1 . The current sphere hits the catch net or a water level before 
encountering any other object. If it does, it is placed there. 

[00151] • 2. The current sphere loses contact with one of its two 
contacts before encountering any other object. If it does, FOBJ2 is recalled. 

[00152] • 3. The current sphere loses contact with both objects 
simultaneously before encountering any other object. If it does, it is in free fall and FOBJ1 is 
recalled. 

[00153] • 4. The current sphere encounters one or more additional 
objects. If it does, STABLE is called to check whether three or more are in compressive contact. 
If they are, the current sphere is placed there. If not, contacts are lost with those objects in 
tensile contact. If two objects remain in tensile contact, FOBJ3 is recalled; if one object, then 
FOBJ2 is recalled; and if all contact is lost, the current sphere is in free fall and FOBJ1 is 
recalled again. 

Detail of Subroutine STABLE 

[00154] The structure and manner of operation of the STABLE subroutine according to 
the preferred implementation will now be addressed in greater detail. Particularly, the STABLE 
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routine may consider the current sphere with weight W acting in the -z direction, if it is touching 
one or more objects, if it is compressively touching them, or if it would roll away from the one or 
more objects. 

[00155] Index numbers of placed pack spheres being touched are passed to subroutine 
STABLE in an array called INDEXTOUCH. A separate parameter called TOUCHWALL states 
whether the cylinder wall of the current sphere's submode is being touched. The manner in 
which the preferred implementation determines the compressive touching criteria when the 
current sphere touches (a) one sphere, (b) two spheres, (c) three spheres, (d) one sphere and the 
wall, or (e) two spheres and the wall will now be addressed. There is generally no need to 
consider compressive touching for more than three objects because the case for more than three 
objects can be broken into sets of three. At least three objects should be checked, however, 
because that is the minimum number needed for stable placement in the presently preferred 
implementation. 

Compressive touching of one sphere 
[00156] Determining compressive touching on one and only one pack sphere can be 
done as follows. If the current sphere is touching its northern hemisphere, it is compressive. 
Because current spheres are dropped from above, they cannot initially touch southern 
hemispheres and FOBJ1 rejects a sphere as a touching object if the current sphere only touches 
its equator. 

Compressive touching of two spheres 
[00157] Consider a current sphere that has descended along a longitudinal line of 
Sphere 1 until it reaches Sphere 2. When it touches Sphere 2, there are only two possibilities: it 
will either roll along a roll corridor formed by compressively touching both Spheres 1 and 2 or it 
will over-roll Sphere 2 and lose contact with Sphere 1. To check for compressive contact with 
both Spheres 1 and 2, the preferred implementation checks whether the current sphere would 
contact Sphere 1 if it were to roll down the longitudinal line of Sphere 2 where the current sphere 
had just contacted Sphere 2. If the current sphere would contact Sphere 1, then it must be in 
compressive contact with Sphere 1. It is certainly in compressive contact with Sphere 2 because 
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of the contact therewith. Therefore, it is in compressive contact with both spheres. If, on the 
other hand, the current sphere would roll away from Sphere 1 while rolling down the 
longitudinal line of Sphere 2, then it will leave Sphere 1 and is only in compressive contact with 
Sphere 2. 

[00158] Consider the positions of the three spheres; (Xj, Yj 9 Zj) for Sphere 1, (X 2 , Y 2 , Z 2 ) 
for Sphere 2, and (X C9 Y c , Z c ) for the current sphere. The radii of these three spheres are aj> a 2 , 
and a c , respectively. The Sphere 2 longitudinal line at which the current sphere is touching 
Sphere 2 is at azimuth 

y c2 = ten-\{Y c -Y 2 )IX c -X 2 )} (16) 

and it touches at latitude 

Q C 2 = cos-X(Z c -Z 2 )/a c -a 2 )]. (17) 

[00159] If 0 c2 > 7i/2, then the current sphere has rolled into the underside or southern 
hemisphere of Sphere 2. In this case, it will be in compressive contact with the two spheres. If 
0 C 2 < n/2, then we must see if it contacts Sphere 1 if it were to roll down the longitudinal line at 

<Pc2. 

[00160] Let the distance between the current sphere and Sphere 1 be denoted D c j which 
is a function of 0 C 2. Its square is given by 

D* = (X c - X x f + (Y c - Y x f + (Z c - ZO 2 = [{a c + a 2 ) sin 9 c2 cos <p c2 +X 2 - X x f + 
[{a c + a 2 )sin 6 c2 sin <p c2 + Y 2 - Y x f + [(a c + a 2 )cos 0 c2 + Z 2 - Z,] 2 . 

If D c i or its square increase as 9 C 2 increases (i.e., as the current sphere rolls down Sphere 2's 
longitudinal line), then it has lost contact with Sphere 1. That is, if its derivative with respect to 
0 C 2 is greater than or equal to zero, it is not in compressive contact with Sphere 1. If the 
derivative is negative, it is in compressive contact with both spheres. The expression below is 
the derivative to within a positive constant of 2(ac + a 2 ), so if it is negative, the current sphere 

37 



compressively touches both Spheres 1 and 2: 



[(a c + a 2 )sin 9 c2 cos <p c2 +X 2 - X x ] cos B c2 cos <p c2 + [(a c + a 2 )sin 0 c2 sin <p c2 + 
y 2 - ^1] cos 6 c2 sin #> c2 - [(a c + # 2 ) cos 0 c2 + Z 2 - ZJ sin 0 c2 . 

[00161] These two-sphere calculations are done in this implementation by Subroutine 



[00162] In the case of compressive touching of three spheres, the current sphere at 
position (X c , Y c , Z c ) and with radius a c is in contact with three pack spheres at positions (Xi, Yj, 
Zj), (X2, Y2, Z2), and (X3, Y3, Z3) and with radii aj, a2, and as. Let the current sphere have weight 
W due to gravity or similar force. Assume it is attached to the three spheres at the point of 
contact with them. The preferred implementation solves a force equation to determine whether 
the attachments are compressive or tensile. If an attachment is tensile, the current sphere will 
roll away from it. 

[00163] Let the unit vectors pointing from the current sphere toward Sphere i where i is 
1, 2, or 3, be denoted by e,. Let the magnitude of the force applied on the current sphere by 
Sphere i be denoted by _/J. If ft is positive, the attachment is tensile; if negative, it is compressive. 
The three equations needed to solve for the three ft are given by the vector equation 



[00164] Because in the implementation we are only interested in whether the fs are 
compressive (fi < 0) or tensile (f t > 0) and are not interested in their magnitudes, the preferred 
implementation simply sets Wio unity and obtains for the / solutions 



ST2S. 



Compressive touching of three spheres 



fie 1 + f 2 e 2 + f 3 e 3 = (0,0-W). 



(20) 



// = {fihfiiy ~ e 2x e 3y )/det E 



f 2 = (e Ix e 3y - e 3x e Jy )/det E 



(21) 
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f 3 = ie 2x e ly - e Jx e 2y )/dQt E 



where det E is the determinant of the 3 x 3 matrix whose columns are the unit vectors eu e 2 , and 

e 3 : 

det E = e lx (e 2y e 3z - e 3y e 2z ) + e 2x (e 3y e Jz - e Iy e 3z ) + e 3x (e Iy e 2z - e 2y e }z ) . (22) 

[00165] If the determinant is zero, the set of equations has no solution because the three 
placed pack spheres and the current sphere all lie in the same plane. In the unlikely event that 
this should happen in practice, the preferred implementation first sees if another set of three 
spheres gives stability if the current sphere is touching more than three placed pack spheres. If 
there is no stable triplet, then it enacts a subroutine to perturb the current sphere's position 
slightly in a downward direction such that it will not overlap any of the three spheres. If the 
perturbation causes contact to be lost with one or two of the touching spheres, it can proceed 
with FOJB2 or FOJB3. If no contact is lost, then at least the four spheres are no longer coplanar 
and it can redo the processing without a zero determinant. 

[00166] These three-sphere calculations are done in the preferred implementation by 
Subroutine ST3S. 

Compressive touching of one sphere and the wall 
[00167] In the case where there is compressive touching of one sphere and the wall, the 
wall cannot be the first object contacted because in this implementation, the current sphere drops 
parallel to it when free falling. Therefore, the only way the current sphere can be in contact with 
both its cylinder wall and one placed pack sphere is if it hits Sphere 1 and subsequently rolls into 
the wall. The current sphere is always in compressive contact when rolling on one sphere. 
Furthermore, because the current sphere must be in contact with Sphere l's northern hemisphere 
when it contacts the wall, it must be in compressive contact with the wall because it just rolled 
into it. Therefore, any time the current sphere is in contact with one sphere and the wall, it is 
always in compressive contact with both and will descend along a roll corridor defined by both. 
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In the unlikely event that the current sphere hits the cylinder wall and nothing else at the same 
time it reaches the equator of Sphere 1, it free falls from that point and FOBJ1 is recalled. 



Compressive touching of two spheres and the wall 



[00168] The case of compressive touching of two spheres and the wall is a simple 
offshoot of the three-sphere case. The cylinder wall will replace one of the spheres as a contact. 
However, in the absence of friction, the wall can only apply a force to the current sphere directed 
radially inward. 

[00169] Let ri and r 2 denote the positions of the two pack spheres and r c denote the 
position of the current sphere. The relative positions of the two pack spheres to the current 
sphere is r ic = r, - r Cy i = 1 or 2. These relative vectors point from the current sphere center to the 
z-th pack sphere center. Let e ic be the unit vector pointing along r lc . The force the current sphere 
applies to each pack sphere and the wall must sum to its weight W. That gives the equation 



[00170] The unit vector e wc , which points radially outward from the cylinder center, is 
given by 



[00171] As before in the three sphere case, the preferred implementation solves for the 
fi, i = 1, 2, or w, but actually only their sign is needed: positive, negative or zero. Hence Wis set 
to unity. Let e ix be the x component of the vector e iC9 i = 1, 2, or w>, and similarly for the;/ and z 
components. Then the solutions for the fi are given by 



Fie }c +f 2 e 2c + f w e wc = 



(0, 0, -KPO- 



(23) 



e wc = (x c \ + rj)/(x c 2 +y c 2 ) 



(24) 



/; = (e 2y e wx - e 2x e^J/dct e 
f2 = (eix ewy - e {y *?wx)/det e 
fw = (e }x e 2y - e\ y e 2x )iteX e 



(25) 
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where 



det e = ~e !x e 2z + e ly e 2z e wx + ei z {e 2x - e 2y e wx ). (26) 

If an/;, i = 1, 2, or w, is negative, the force between the current sphere and that object is 
compressive and they stay in contact. If it is positive, the force is tensile and the current sphere 
will roll away from that object. In the unlikely event that it is zero, the implementation considers 
it tensile and has the current sphere roll away from the object, thereby trying to find a point of 
stability at lower altitude. If all three are negative, the sphere is in stable contact and is placed 
at its current position. 

[00172] These stability calculations for two spheres and a wall are performed in the 
preferred implementation by Subroutine ST2S1W. 

Detail of Subroutine FOBJ1 

[00173] Returning to the subroutine FOBJ1, it finds the first object the descending 
current sphere encounters as it is "dropped" from the position selected by subroutine 
DROPSITE. There are three possible objects the current or selected sphere can encounter: (1) an 
already-placed pack sphere, (2) a water level, or (3) the catch net. A fourth object that may be 
encountered later is the cylinder wall, but because the current sphere is assumed to drop parallel 
to the cylinder wall, it cannot contact the wall until it rolls off some placed pack sphere. 

The catch net 

[00174] In accordance with the method aspect of the invention, the selected particle 
placement includes defining a catch net representative of buoyancy of a portion of the placed 
particles and positioning the catch net within the space based upon the placement of the portion 
of the placed particles. The selected particle is placed in non-overlapping relation with respect to 
the catch net. Preferred versions of the method include defining a pack surface for the placed 
particle, and positioning the catch net relative to the pack surface. The pack surface may 
comprise a top layer of the placed particles within the space, or preferably, within each of the 
subspaces. The particle surface optionally may correspond to an average of the south poles or 
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south pole positions of the top layer of placed particles. The catch net may be positioned at a 
distance from the pack surface based upon a selected particle radius, such as an average particle 
radius of the top layer or another subset of the placed particles at the pack surface. 

[00175] In the way of background, when simulating a particle pack, particularly one 
with a wide particle size distribution, if the volume fraction of large particles is much larger than 
that of small particles, the small particles can rattle down through the pack by falling through the 
interstitial spaces of the large particles. In many particulate materials, there are buoyancy or 
suspension considerations that would keep the small particles more or less well distributed 
throughout the pack. Examples of such buoyancy or suspension factors may include the nature 
or properties of the binder or other inter-particle medium or media, the interaction of the particles 
with the inter-particle medium or media, interactions among the particles themselves (such as 
electrical or magnetic interactions), interactions with external forces or fields that tend to 
disperse or suspend the particles, and similar phenomena. 

[00176] In accordance with this aspect of the present invention, this buoyancy or 
suspension phenomena is simulated using a "catch net." The catch net comprises an assumed 
surface or construct to catch the current or selected particle as it is placed and to prevent that 
particle from passing a defined location or level which theoretically represents space occupied by 
smaller ones of the previously placed particles. 

[00177] In the presently preferred implementation, a catch net is incorporated into the 
space, below which the south pole of the falling particles cannot descend. When a particle's 
south pole hits the catch net belonging to its cylinder, it is placed in the pack at that point, even if 
it is not touching any other particle. Preferably there is a separate catch net for each particle 
submode, and thus for each subspace or cylinder. Each catch net preferably covers the entire 
cross section of that submode's cylinder. If the current sphere's south pole hits that catch net, it 
is placed there. It does not see or respond to any other submode's catch net. The altitude of the 
catch net preferably is determined dynamically, that is, it can depend on which placed particles 
are being contacted. 

[00178] The use of a catch net is not necessarily limited to a specific type or size of 
particle, or a particular size or geometry of the space or subspaces. In accordance with the 
presently preferred implementation, however, as noted, the particles comprise spheres and the 
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characteristic dimension of each of the particles comprises a radius. The central string comprises 
a line, preferably the z axis, within the space. The space comprises a cylinder, and the subspaces 
comprise concentric cylinders disposed about the central string, each of the cylinders defining 
one of the subspaces and having a radius that corresponds to the radius of the particles in the 
category corresponding to that cylinder. 

[00179] In accordance with the presently preferred implementation, the method further 
includes defining a pack surface for the placed particles, and the catch net positioning comprises 
positioning the catch net relative to the pack surface. Also in accordance with this 
implementation, the catch net positioning comprises positioning the catch net a distance from the 
pack surface based upon a selected particle radius. 

[00180] The placed particles may be assumed to comprise a top layer, and the pack 
surface optionally but preferably comprises an average of the particle locations of the placed 
particles in the top layer. In this implementation, the particle surface corresponds to the south 
poles of the top layer placed particles and more specifically corresponds to an average of the 
south pole positions of the top layer placed particles. 

[00181] In the presently preferred implementation, the catch net extends across the cross 
sectional area of the space, across the cross sectional areas of each of the subspaces. 

[00182] In this preferred implementation, the catch net comprises N subnets, one of the 
subnets corresponding to each of the subspaces. Each of the concentric cylinders 1, 2, 3, ... N 
has a circular cross section perpendicular to the central string. In the preferred implementation, 
the catch net comprises a plurality N of subnets, or individual catch nets, one for each of the 
cylinders. Each of the subnets preferably extends over the cross sectional area of the 
corresponding subspace or cylinder. The position or level of the catch net or subnet for each of 
these subspaces can and typically will differ from one another, although in some instances they 
may be at the same or nearly the same level. 

[00183] Preferably, the positioning of the catch net comprises setting the catch net or 
each of the subnets at a selected distance from the top surface. This selected distance may be a 
user-defined quantity, or it may be dictated by the material being simulated or other factors. The 
catch net preferably is positioned at an altitude along the z-axis that is a function of the particle 
radius of particles corresponding to the mode or submode of the cylinder in which the catch net 



43 



or subnet resides. The catch net or subnet preferably, but optionally, is some multiple of the 
particle radius, preferably one or two ratio, below the top layer or pack surface. Adjustments or 
offsets, however, may be used. 

[00184] When the pack is first being formed so that not enough particles have yet been 
placed to adequately define a pack surface, the catch net may be defined as the bottom of the 
cylinder into which the particles are being dropped. In some instances, however, particularly 
where the specific application involves simulating a very large three-dimensional pack, it is often 
advantageous or desirable to avoid the edge effects of the flat bottom of the cylinder. 
Accordingly, in some instances the positioning of the catch net may comprise spacing the catch 
net from the base surface of the space. The spacing from the base surface, for example, may 
simulate the positioning of the catch net for a top layer of the placed particles, for example, the 
positioning of the catch net off of the base of the space may be done in a way that the bottom 
layer of placed particles simulate the shape, configuration and appearance of the top layer of 
particles after the particle pack has been constructed to include multiple layers. 

[00185] This may be implemented by assigning as an initial catch net position for a &-th 
one of the subspaces Zj n i t (k), assigning as the characteristic dimension of the particles of a k-th 
one of the particle categories ak, assigning as the characteristic dimension of a small one of the 
particles a m j n of the particles and assigning as the characteristic dimension of a large one of the 
particles a max . Below a threshold level the catch net positioning further comprises positioning the 
catch net for a k-th one of the particle categories at a catch net position Zn e t(k) within a Ar-th one 
of the subspaces determined by 

Znet(k) = Zi n u + H r a k a min /a max (27) 

where r represents a weighting coefficient and H represents a switching coefficient. The 
weighting coefficient preferably, but optionally, is assigned a random number. Below the 
threshold value the switching coefficient preferably is assigned a value of one, and above the 
threshold value the switching coefficient preferably is assigned a value of zero, 

[00186] In the presently preferred implementation, the switching coefficient mechanism 
is implemented by a switch called ROUGHNET which, when true, invokes the random number 
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generator to slightly change the altitude of the catch net so that the first set or layer of placed 
particles do not all have the same altitude. If ROUGHNET is false, the south poles of the first 
set of placed particles do have the same altitude because they hit the same catch net. 

[00187] Let Zi n u be the initial altitude above which the pack will be built. It is specified 
by the user in the preferred implementation but its default preferably is zero. The initial catch 
net altitude for the &-th mode is 

Znet(k) = Z init + H r a min /a max (28) 

where r is a random number in [0, 1] a min is the smallest particle radius of all submodes 
belonging to the pack, a max is the largest, and H is unity if the ROUGHNET switch is true and 
zero if it is false. This expression for the catch net altitude holds until enough particles have 
been placed to establish a pack surface, i.e., until m (which is greater than or equal to m min ) is 
reached. No part of any particle is permitted to descend below Z init in this implementation. 

[00188] The catch net preferably is positioned relative to the pack surface. The pack 
surface can be determined in a number of ways. It may be determined, for example, based on the 
large particles in the top layer, or another subset of the placed particles, e.g., such as the last 
particles placed. 

[00189] The particle surface may be ascertained in the following manner. For each 
particle category k of the placed particles in the top layer, the program defines a particle radius aj 
for the placed particles i of that category k. For the cylinder k corresponding to the particle 
category k, the program assigns a cylinder radius W k . The program also assigns a top layer 
particle number m(k) and determines values for m(k) by evaluating 

m{k)-\ m(k) 

I a* <W k 2 <Zal ,k = 1,2,.-, N 

i=\ i=\ y^y) 

Sub mod e(i)<Jc Sub mod e(i)Zk 

where N is the number of particle categories, e.g., the number of submodes. The program then 
determines the pack surface location using 
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1 m 

S = -Z(Zi-ad 



(30) 



where S represents the pack surface and Z\ represents the position of a center of a center one of 
the placed spheres. 

[00190] In accordance with the preferred implementation, the surface altitude S of the 
pack is defined as the average south pole altitude of the last m placed particles where m is a 
positive integer that will be discussed below. This surface altitude is given by 



where Z, is the altitude of the center of Sphere i which is one of the last m placed particles and a t 
is its radius. Instead of using a single i-th particle radius in Eq. 31, one may generalize the 
expression by permitting {Z x - aO to be written {Z\ - f(aj)) where f(aO is a function of particle 
radius aj. The preferred value for f(aO = aj, but f(aj) may assume other forms, for example, such 
as f(ai) = nai where n is a positive real or integer number. 

[00191] Although the "surface" of a random particle pack may be defined in a number 
of ways, the surface preferably is defined only after the integer m is large enough that the falling 
spheres have covered the cylinder area. The A:-th cylinder with radius W k has area n W k 2 . Let m 
be defined as the number of most recently placed spheres whose submodes have cylinder radii 
less than or equal to W k and whose summed cross-sectional area just exceeds or equals the &-th 
cylinder area. For the &-th submode, m is taken from the following condition. Let m(k) be the 
number m for the &-th submode. It is found from 



1 m 

S = -Z(Z i -a l ) 



(31) 



m(k)-l m(k) 



E a] <W k 2 < I aj ,k = \,2,...,N 

i=\ i=\ 



(32) 



where N is the total number of submodes. The submode k for which this condition is satisfied 
first, that is, has the lowest m(k) becomes the m for the pack: 
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m = MtN[m(k)]. 



(33) 



[00192] In the preferred implementation there is a minimum number m min for those 
unused packs that defeat the reasoning leading to the value for m given in Eqs. 32 and 33. The 
default value of m m i n may be assigned to a value of 12 in the preferred implementation. 

[00193] Once a surface is established, its altitude S is determined by using Eq. 31. As 
the current sphere descends, if it hits a pack sphere much larger than itself, and if the TUNNEL 
switch is true, it will tunnel through the large sphere to check for possible cavities under the large 
particle. Therefore the catch net should be placed quite low so the current sphere can investigate 
possible cavities. If TUNNEL is false, the catch net does not need to be that low. 

[00194] If the current sphere hits a sphere somewhat bigger than itself, then it is unlikely 
to tunnel because the surface sphere is not large enough to accommodate it in a cavity. Then the 
catch net need not be as deep. 

[00195] Finally, if the current sphere hits a sphere smaller than itself, then it cannot 
penetrate too deeply into the pack so the catch net can be rather shallow. 

[00196] Letting the subscript i denote the first pack sphere contacted, the catch net 
criteria are quantitatively summarized as follows: 

[00197] If ai/a c < 1, then the catch net is shallow: 



ZnetQt) ~ S - <Xi. 



(34a) 



[00198] If 1 < aila c < a x , where a x preferably is V6 + 2, then the catch net is moderate: 



Z net (k) = S-2a c . 



(34b) 



[00199] If a/a c > a x , where a x preferably is V6 + 2, then the catch net is deep: 



Znet(k) = -2a c - at. 



(34c) 
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[00200] The value of a x represents the sizes required for a current sphere to fit into a 
cavity formed by larger spheres, which preferably assumes a value of (V6 + 2). If the centers of 
four identical larger spheres form a regular tetrahedron, and a smaller sphere just fits into the 
interstitial space of the tetrahedron, then the ratio of the smaller sphere to larger sphere diameter 
is V6 + 2 in preferred implementation. 

[00201] Further in accordance with the preferred implementation, the code monitors the 
selected or current particle placement to identify the occurrence of the following sequence of 
events: 

[00202] 1. The catch net is set at an altitude when the current sphere hits the 

z-th pack sphere. 

[00203] 2. Then the current sphere rolls off the z-th sphere and loses contact 

with it. 

[00204] 3. When it hits another pack sphere, the catch net is set above the 
current altitude of the current sphere's south pole. 

[00205] If this series of events is detected, the current sphere is raised (its position is 
adjusted upward) until its south pole is at the new catch net altitude, and it is then placed there. 

The water level 

[00206] In accordance with another aspect of the invention, the placement of the 
particles comprises defining a water level representative of a level of a portion of the placed 
particles that are smaller than the selected particle and represent a surface of the smaller placed 
particles, and positioning the water level within the space based upon the smaller particle 
surface. The current or selected particle is then placed in non-overlapping relation with respect 
to the water level. 

[00207] The water level serves a different purpose than the catch net. Recall that the 
particle pack according to the preferred implementation is constructed using a series of 
concentric cylinders, one cylinder defined for each particle category or submode. Recall further 
that the centers of particles in this implementation cannot pass outside their own corresponding 
and constraining cylinder into which they are dropped. Large cylinders are used for large 
particles and small cylinders for small particles, unless the user for some reason designates 
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otherwise. The cylinder construct facilitates or permits use of the reduced dimension approach 
to particle pack construction. Its purpose is to prevent the program code from having to fill large 
cylinder volumes with potentially very large numbers of small particles. 

[00208] This approach, however, creates the problem that the large cylinders will 
contain empty space into which large spheres can roll during simulated particle placement where, 
in a corresponding physical particle pack such as the one being simulated, this would not occur. 
In the actual pack, the smaller particles would be included even in the regions comprising the 
larger cylinders, and these smaller particles would be encountered to prevent such rolling. 

[00209] The water level simulates this presence in the larger cylinder areas of smaller 
particles that, according to the presently preferred implementation, have not been placed in the 
larger cylinders. 

[00210] There are a number of methods or approaches for selecting or setting the water 
level. According to one such approach, the water level positioning comprises determining an 
average location of the particle locations along the central string of the particles of the portion of 
placed particles that are smaller than the current or selected particle, and positioning the water 
level at the average location. The water level may be set at the average altitude of the small 
particle "surface" as if the outside cylinders were filled with small particles. Then when a large 
particle or sphere descends onto the water level, the particle, and preferably its south pole, is 
stopped there as it would be if that volume had been filled with the smaller particles. 

[00211] In the presently preferred implementation, wherein the plurality of sub spaces 
(here concentric cylinders) is used, the water level comprises a plurality of subspace water levels, 
preferably one subspace water level for each subspace or cylinder. These water level positions 
may and typically will differ from one another as the particle pack is constructed. The water 
level positioning optionally, but preferably, comprises assigning a subspace water level position 
to a portion of the subspace water levels. The portion of the subspace water levels preferably 
correspond to one water level for each cylinder except the first cylinder. In the preferred 
implementation, all of the cylinders except the first or central one, containing the smallest 
particles, has a subspace water level and corresponding subspace water level position. 

[00212] In the presently preferred implementations, the water level positioning 
comprises using an offset to position the water level, or an offset relative to the average location 
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of the particle locations. Preferably, the water level positioning comprises using an offset for 
each of the subspace water level positions. The water level positioning similarly may comprise 
using an offset for each of the subspace water level positions. 

[00213] In the presently preferred implementation, the water level altitude for the i-th 
cylinder (i> 1) is defined as the average exposed north pole altitude of all particles with radius 
less than ai. The user may raise or lower this level by adding a user-defined offset, but the 
default for this offset in this implementation is zero. This water level exists between the i-th 
cylinder wall of radius Wi and the next smaller cylinder wall of radius Wi-j. Thus the spaces 
between different cylinder walls have different water levels. The cylinder corresponding to the 
smallest particle's submode has no water level, only a catch net. 

[00214] If a larger sphere is descending onto the pack near the center cylinder, it may hit 
the center cylinder's pack of small particles and tend to roil off its edge and fall to the floor of 
the pack or fall until it hits another large sphere. Thus the pack would segregate with larger 
particles outside smaller particles. With the water level in place, however, if the larger sphere 
rolls off the pack in the central cylinder, it simply "floats" on the water level which is at the same 
level as the "surface" of the central cylinder pack. Thus, any time a descending sphere hits a 
water level, it is placed where it is, just as was done with the catch net. Ordinarily, water levels 
will be higher than catch nets, although this is not always true. 

Search for the first object encountered 

[00215] A free- falling current sphere can encounter the following situations. It might hit 
a single pack sphere. If so, FOJB2 is called. It can hit two pack spheres simultaneously. If so, 
FOBJ3 is called. It can hit three or more pack spheres simultaneously. If so STABLE is called 
to see if it is in stable contact with any triplet of these spheres. Hitting two or more pack spheres 
simultaneously is highly improbable in a double precision random code but not so improbable if 
the pack is semi-ordered, and hence these scenarios are included in the preferred implementation. 

[00216] In addition, the south pole of the current sphere may hit the catch net or a water 
level before it hits any other sphere. In either case, it is placed where it is. 

[00217] A pack sphere, e.g., the y'-th sphere, is a candidate for being the first object 
encountered by the descending current sphere if their two radii overlap, i.e., if 
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(Xc-Xj) 2 + (Y c -Yj) 2 <(a c + aj) 2 



(35) 



If there were no other obstacles in the way, the current sphere would strike Sphere j at an altitude 
of 

Zc{j) = Zj + [(a c + ajf - (X c -Xj) 2 - (Y c - Yjff 2 . (36) 

[00218] Once the pack has grown to many spheres, efficiency can be greatly enhanced 
to the extent that as many spheres as possible be excluded from the detailed considerations of the 
above two equations, while not unduly adversely affecting the accuracy of the simulation in 
representing the physical particle pack being simulated. For example, if | x c - xj | or | y c - yj \ 
alone were greater than (a c + aj), then there is no need to check the square of their sums as 
above. 

[00219] The processing demands are reduced in the preferred implementation because 
all spheres whose north poles have submerged beneath the catch net are off limits and need not 
be included. Simply flagging them by putting a minus sign in front of their mode or submode 
number indicates that no calculations are to be done to see if they are in range of the descending 
current sphere. Furthermore, in this implementation the program keeps track of the lowest index 
number for any sphere whose north pole protrudes above that mode's catch net so that the 
eventually large list of particles beneath the catch net are not searched at all. 

[00220] Submergence of a placed pack sphere beneath the water level, however, is not 
necessarily grounds for dismissing it from consideration. The current sphere may roll into a 
cylinder for which the water level is different and the pack sphere may again become a candidate 
for contact. 

[00221] If the current sphere is dropped far enough outside a set of cylinders that there is 
no chance for a first contact with the spheres belonging to those cylinders, then in the preferred 
implementation such sphere(s) are dismissed from further consideration. If the submode number 
of a sphere belongs to an out-of-range cylinder, no further processing is required in the preferred 
implementation. Let Wk be the cylinder wall radius for the k-th particle submode. If 
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(W k + a k ) 2 <X c 2 + Y c 2 



(37) 



then no k-th mode sphere is a candidate for first contact with the descending current sphere in 
this implementation. 

[00222] In summary then, FOBJ1 searches the list of all placed pack spheres above the 
catch net for the submode to which the current sphere belongs, excluding all placed spheres that 
may be in interior cylinders that the current sphere cannot contact or overlap. The possibility of 
overlap in the x directions is checked, followed by the possibility of overlap in the y directions, 
followed by a direct calculation of overlap. The altitude z c (j) at which the current sphere would 
contact Sphere j is stored and that sphere is rejected if other placed pack spheres have a higher 
contact altitude. Again, if no pack spheres make contact before the current sphere reaches the 
catch net or water level, it stops and is placed there. 

Detail of Subroutine FOBJ2 

[00223] Subroutine FOBJ2 finds the second object to be contacted by the descending 
current sphere. FOBJ2 would only have been called if the current sphere had first contacted a 
single pack sphere that has been denoted Sphere 1. If it had contacted the catch net or water 
level, the current sphere would have been placed there and FOBJ2 would never have been called. 
If it had contacted more than one sphere simultaneously, FOBJ3 and/or STABLE would have 
been called. FOBJ2 is called only if a single pack sphere's northern hemisphere has been 
touched. 

[00224] When the current sphere touches the northern hemisphere of Sphere 1, it begins 
to roll down the path of steepest descent which is Sphere l's longitudinal line passing through the 
point of contact. The volume swept out as the current sphere descends along this path toward the 
equator is called the roll corridor. The end of the roll corridor is either where the current 
sphere's center reaches the equatorial plane of the Sphere 1, reaches its cylinder wall, or reaches 
the catch net or water level, whichever is higher. The preferred implementation identifies all 
other pack spheres that protrude into the roll corridor. 
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[00225] Processing can be somewhat simplified if the cylinder wall is assumed to be the 
boundary that the center of the current sphere cannot pass, rather than the boundary through 
which its leading edge equator cannot pass. 

[00226] If the center of the current sphere intersects the cylinder wall during descent, 
then that point is the end of the roll corridor because it will at least have found the wall as its 
second object if no other pack sphere has been contacted first. 

[00227] The beginning of the roll corridor is then the place of initial contact between the 
current sphere and Sphere 1. The current sphere's polar angle at that initial contact is given by 

9 ci = cos\(Z c - Zj)/(a c + aj)]. (38) 

[00228] The end of the roll corridor 0 cf will either be at the equatorial plane of Sphere 1 
or where the center of the current sphere intersects the catch net, the water level, or the cylinder 
wall. If it is at the equatorial plane, then 

9 cf = 7c/2. (39) 
If the altitude of the catch net Z net >Zj- a c , then 

e C f = cos -1 [(Znet + a c - Zj)/(a c + ai)] . (40) 
If the altitude of the water level Z w >Zj - a c , then 

9cf = cos J [(Z w + a c - Zj)/(a c + ai)]. (41) 

[00229] The last possibility is that the roll corridor ends with an intersection of the 
cylinder wall. Let <p c be the azimuthal position at which the current sphere has touched Sphere 1. 
It is given by 



9c = tan- , [(y c -y y )/(Z c -Jr / )]. (42) 
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[00230] The radial position of the current sphere as it descends Sphere 1 as a function of 
9 C is given by 



{[Xi + (a c + ctj) sin0 c coscp c ] 2 + [Yi + (a c + aj) sin9 c sin9 c ] 2 } 



(43) 



[00231] If this radial position is equal to the cylinder radius W for some value of 9 C 
between 0 and nil, then the roll corridor intersects the cylinder wall. Setting this expression 
equal to W, squaring both sides, and solving the quadratic equation for sin9 c gives 



[00232] Because 9 cf is a polar angle, it lies in the range [0, n]. Therefore, sin9 c f can 
never be negative. If both right-hand side solutions of Eq. 44 are greater than unity, there is no 
solution and the current sphere does not intersect the cylinder wall. If both are less than or equal 
to unity, then we choose the larger one. 

[00233] The presently preferred implementation conducts appropriate processing to 
evaluate each of these possible circumstances. Before evaluating Eq. 44, the program code 
according to this implementation adds the global radial position of Sphere 1 to the sum of the 
radii of Sphere 1 and the current sphere. If this sum is less than or equal to W, then the sphere 
cannot hit the cylinder wall in any case so the more involved evaluation above is unnecessary 
and preferably is not performed. 

[00234] The end of the roll corridor is the minimum 9 cf of Eqs. 39, 40, 41, and 44. If the 
current sphere encounters the catch net or water level, it is placed there. 

[00235] The preferred implementation also checks for all pack spheres whose volumes 
intrude upon the roll corridor swept out by the current sphere which covers the arc from 9 ci to 9 cf 
and azimuth cp c on Sphere 1 . If there is a 9 C with smaller value than the end of the roll corridor 
9 c f at which the current sphere encounters one or more pack spheres, then the sphere or spheres 



sin 0 C/ = 



-(X, cos <p c + Y ] sin <p c ) ± ^J(X X cos (p c + Y x sin <p c ) 2 + W 2 - X x 2 - Y x 2 

a c + a x 



(44) 
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touched are noted. Program control is then passed to subroutines that deal with simultaneous 
contact of two or more spheres, as described herein. 

[00236] If no pack spheres overlap the roll corridor, then the current sphere is placed if it 
hits the catch net or water level. If it hits neither of these, then it rolls to the equator and drops 
from there. At that point it is a free falling sphere again so FOBJ1 is recalled. 

[00237] We now turn our attention to finding which placed pack spheres overlap the roll 
corridor and at what point the current sphere would strike them as it descends from 6 C j. The term 
"global coordinate system" is used herein when discussing a single coordinate system in which 
the entire pack is defined. Capital letters are used when referring to a global coordinate system, 

e.g., R; = (Xj, Y t , Zi) for the Cartesian position of the /-th sphere. 

[00238] A local coordinate system, on the other hand, is one whose origin is at the center 
of one of the spheres and whose axes are parallel to the global axes. We use small letters to 
define it, e.g., r t = (*,-, y u z,) for a position relative to the center of the /-th sphere. 

[00239] Consider the y-th pack sphere. Its center is at global position Rj = {Xj, Yj, Zj) and 
its radius is aj. In the global coordinate system, the position of the current sphere as it descends 
along Sphere / is 

R c - Ri + (a c + a,) (sin 9 C cos q> c i + sin 0 C sin <p j + cos 0$.) . (45) 

[00240] The current sphere rolls into Sphere j if there is a solution 9 C that lies in the 
interval 0 c j <0 C < 0 c f to the equation 



R c - Rj =a c + aj . (46) 



[00241] If this equation can be satisfied at all, it will have two solutions: one where the 
current sphere touches Sphere j from above and another where the current sphere touches Sphere 
j from below. (If it just nicks Sphere j, these two solutions can be at the same 6 C , but if Sphere j 
is just nicked, the preferred implementation considers that a miss.) We are only interested in the 
top or smaller 6 C solution, and then only if the top 6 C lies in the roll corridor range. 
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[00242] Explicitly, the equation to be solved is 

[Xi -Xj + (a c + ai) sin0 c coscp c ] 2 + [Y t - Yj + (a c + a t ) sin9 c sincp c ] 2 
[Z/ - Zj + (a c + A,) cos9 c ] 2 = (a c + a,-) 2 

which, defining Xy = Xi- Xj 9 Yy = 7/ - )y, and Zy = Z; - Z,, can be rewritten 

(A!lj cos(p c + Yy sincpc) sin0 c + Zy cos0 c = [(fl/ - a){aj + a,- + 2a c ) - Ry 2 ]/[2(a c 

where 

i?,/ = Xif + + Zy 2 . 

This equation has the form 

A sin9 c + 5 cos0 c = C 

where 

^ = Xy coscpc + 7y sincpc 

B = Zy 

C = [(a, - fl/Xfl,- + a, + 2a c ) - Ry 2 ]/[2(a c + a,)]. 
Dividing through by V(/i 2 + B 2 ), we obtain 

cosa sin0 c + sina cos0 c = CH(A 2 + 5 2 ) 

where 

cosa-^(/V(^ 2 + 5 2 ) 
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and 



sim=BN(A 2 + B 2 ). (56) 
[00243] A common trigonometric identity turns Eq. 54 into 

sin (0 C + a) = CN(A 2 + B 2 ) (57) 
which has solution 

6 c = sm A [CH(A 2 + B 2 )]-a (58) 

where 

a = \.zri\AIB]. (59) 

If | CH(A 2 + B 2 ) | > 1, then there is no solution and the j-th sphere does not overlap the roll 
corridor. 

[00244] The FORTRAN function ATAN2 (A 9 B) or its equivalent in other programming 
languages uses information on the signs of A and B to determine in which quadrant a lies so it is 
uniquely defined. The arcsine, however, has two solutions. One solution is the principal value 
which lies in [-nil, n/2] and is found by the standard FORTRAN function ASIN or its 
equivalent. The other solution is 7i minus this principal value. These two solutions comprise the 
two solutions for the current sphere contacting Sphere j on either side along the longitudinal line 
defined by azimuthal angle (p c . Denoting the principal value by capitalizing the arcsine, the two 
solutions are 

e c i = Sin 1 [CN(A 2 + B 2 )] - a (60a) 
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and 



0 c2 = 7t - Sin 1 [CN(A 2 + B 2 )} - a. (60b) 

[00245] If either of these 9 lie outside [0, 2%\, then In must be added or subtracted 
enough times to bring the solutions back into [0, 2%\ If neither of these solutions lie in the open 
interval (0 C i 5 0 C f), then the current sphere does not contact Sphere j. If one solution does lie in 
(9ci, 9cf), that is the position at which the second contact is made. If both do, the lesser of the two 
solutions is the position of second contact. 

[00246] Finally, whichever pack sphere has the smallest 0 C is Sphere 2. That point of 
contact is denoted 0 CJ . If two or more pack spheres are contacted simultaneously, then all their 
index numbers are returned by FOBJ2. If the cylinder wall and one or more pack spheres are 
contacted simultaneously, their index numbers are returned and the flag indicating a wall contact 
is activated. 

[00247] As the preferred implementation searches through the pack spheres to find 
which is touched first, there are several efficiency criteria that can be used to eliminate most 
spheres from the detailed consideration of the above equations. 

[00248] First, if the roll corridor is far enough away from the global Z axis, then the 
current sphere may not be able to approach some inner cylinders closely enough to touch their 
spheres. Therefore, those pack spheres whose mode number indicates they belong to the inner 
cylinders need not be considered further. The preferred implementation can determine which 
cylinder's spheres can never intersect the current sphere's roll corridor. 

[00249] Because the current sphere rolls down a longitudinal line on Sphere z, when 
viewed from above, its projection is a straight line segment on the XY plane with beginning point 
(X C i, Y a ) and ending point {X& Y c j) where 



X ci = Xj + (a c + ai) sin0 ci coscp c (6 1 a) 

Y C i = Yi + (a c + ai) sin0 ci coscp c (61b) 

X cf = Xi + (a c + ai) sin0 cf cosq> c (6 1 c) 

Y cf = Yi + (a c + ai) sin0 cf cos(p c . (6 1 d) 
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[00250] To find whether this line segment comes close enough to the /-th cylinder to 
touch any of its spheres, the preferred implementation turns the current sphere's position into a 
pair of parametric equations: 



X c — X c i + (X c f - X C i)s 
Y c =Y ci + (Y cf -Y ci )s 



(62a) 
(62b) 



where s is a dimensionless parameter that is zero at the beginning of the roll corridor and unity at 
the end. In order for the current sphere to have any chance of touching an /-th mode sphere, it 
must come closer than the sum of its radius plus an /-th sphere's radius from the /-th cylinder 
wall which has radius Wi. In other words, some portion of the roll corridor must come within 
Wi + ai +a c of the global Z axis. This condition is expressed by 



for a value s somewhere in the open interval [0, 1]. If Eq. 63 is never satisfied, then the current 
sphere cannot touch any of the /-th mode spheres. Explicitly written, Eq. 63 is a quadratic 
equation in s: 



(63) 



As 2 + 2Bs + C = 0 



(64) 



where 



A = (X c/ -X ci ) 2 + (Y cf -Y ci ) 2 
B =X cl {X c f- X a ) + Y ci {Y c f- Y CI ) 
C = X cl 2 + Y ci 2 -(W, + a l + a c ) 2 . 



(65) 
(66) 
(67) 



[00251] The solutions are 
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s = 



-B±yjB 2 -AC 



(68) 



A 



If B 2 -AC is negative, there is no solution. If it is zero, there is only one real solution which 
means the roll corridor just grazes the /-th cylinder, in which case it is considered a miss. If it is 
positive, there are two real solutions. The preferred implementation solves for both of them, and 
checks to determine whether either one lies in the open interval [0, 1]. If one or both lie in this 
interval, the /-th mode spheres are candidates for touching the descending current sphere. 

[00252] To find which cylinders are inaccessible to the current sphere, the preferred 
implementation starts with the innermost cylinder (7 = 1) and continues outwardly until it 
reaches a cylinder that is big enough that it is accessible. All larger cylinders also will be 
accessible. 

[00253] There also is a simple altitude criteria to delimit detailed consideration of pack 
spheres as candidates for touching the current sphere. The global altitudes of the beginning and 
end of the roll corridor are 



respectively. Therefore, if a y-th pack sphere lies below the end of the roll corridor by the 
amount 



Z ci = Zf + (a c + a t ) cos0 C i 
Z c /= Zi + (a c + ai) cos0 cf , 



(69a) 
(69b) 



Zj<Z c f- a c - cij 



(70a) 



or lies above the top of the roll corridor by the amount 



Zj > Z ci + a c + a.j 



(70b) 



then it is not a candidate for touching the current sphere. 
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[00254] Likewise, spheres lying beyond the roll corridor in the x and y directions need 
not be further considered, i.e., if 



Xj<X C f- a c - oij 
Yj< Y c f-a c -aj 



(70c) 
(70d) 



or 



Xj>X C i + a c + aj 
Yj>Y ci + a c + a j9 



(70e) 
(70f) 



then the 7-th sphere lies beyond the roll corridor in the horizontal plane and is not a candidate for 
touching the current sphere. 

Detail of Subroutine FOBJ3 

[00255] Subroutine FOBJ3 finds the third object to be contacted by the descending 
current sphere when it is already in contact with two and only two other objects. There are only 
two possibilities for the two initial contacts: (1) both are placed pack spheres or (2) one is a 
placed pack sphere and the other is the cylinder wall for the submode of the current sphere. Each 
of these possibilities will now be addressed in turn. 



[00256] It is presumed at the outset of this discussion that the two touching spheres, / 
and j 9 are in stable (compressive) contact, subroutine STABLE already having been called to 
remove them if they are not. A roll corridor follows a circle defined where the current sphere is 
in contact with both Spheres i and j. Of course, the current sphere cannot touch both spheres 
unless the condition r,, < a,- + aj + 2a c is satisfied where r y7 is the center-to-center distance 
between Spheres i and j. 

[00257] The current sphere's descent along this roll corridor will be defined by a single 
variable cp c '. This variable is the azimuthal angle of a local tilted coordinate system whose origin 



Touching two spheres 
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is in the center of Sphere z, with z axis tilted so it passes through the center of Sphere j, and with 
x axis swung upward so that it lies in a vertical plane. In that coordinate system (which is 
denoted herein by primes), (p c ' alone describes the current sphere's position along the roll 
corridor. In the unlikely event that the current sphere finds itself perfectly balanced on Spheres i 
and j so that q> c ' = 0, then a random number generator determines in which direction the current 
sphere will roll - toward negative cp c ' or toward positive <p c '. 

[00258] To find the coordinates of a pack sphere, say the A:-th sphere, in this primed 
coordinate system, the preferred implementation first translates the global sphere positions to the 
unprimed local coordinate system whose origin is at the center of Sphere i but whose axes are 
parallel to the global axes: 

r ki = R k -R h (71) 

[00259] Capital letters are used herein to denote coordinates in the global coordinate 
system, small letters in the unprimed local system (local to Sphere /), and primed small letters in 
the local tilted coordinate system. 

[00260] The Euler matrix used to rotate the unprimed into the primed coordinate system 
is given by 



A = 



^ - cos Oji cos <pji - cos Qji sin <p }i sin 0 M ^ 

sin <pji - cos (p }i 0 

^ sin cos sin 0,7 sin zo%6 yiJ 



(72) 



where Op and (p^ are the polar and azimuthal angles of Sphere j relative to Sphere i and are given 
by 

cos dji=(Zj-Zi)/rj, (73a) 



and 
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yj^tm'KYj-Yd/iXj-Xd] 



(73b) 



This Euler matrix allows the preferred implementation to obtain the k-th sphere's position in the 
primed coordinate system as 

4* = A. r ki . (74) 

[00261] The current sphere's initial position in the primed coordinate system is given in 
like manner: 

r ci = R c -Ri (75) 
r ci ' = A*r ci . (76) 

[00262] Now that the three components of the vector r ci * have been obtained, the 
preferred implementation can find the initial azimuthal position 

q> ei ' = tm- l (y ei '/x ei '). (77) 

If (p C j' is positive, the current sphere will roll toward higher positive values. If it is negative, it 
will roll toward more negative values (higher in magnitude). In the unlikely case that it is zero, a 
random number generator chooses which direction it will roll. 

[00263] The polar angle of the center of the current spheres in the primed coordinate 
system is constant as long as it is in the roll corridor. It is obtained from the law of cosines as 

cos0 c ' = [r/ + (a c + ad 2 - (a c + aj) 2 ][2 r J( (a c + ai )]' 1 (78) 



with 



sine c ' = +[l -cos 2 9 c '] 1/2 . (79) 
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[00264] There are several possible outcomes as the current sphere descends its roll 
corridor. It can 

• hit the catch net or water level; 

• lose contact with Sphere i or Sphere j or both simultaneously; 

• hit the cylinder wall; or 

• hit another pack sphere or several pack spheres simultaneously. 

[00265] The roll corridor ends when the current sphere encounters the first of these 
situations. The preferred implementation then finds the cp c ' at which one of these events happens. 
The event with the smallest magnitude for <p c ' will be the bottom of the roll corridor and will be 
denoted <p C f'. 

Hitting the catch net or water level 
[00266] First, the preferred implementation finds the global altitude of the current sphere 
as a function of cp c ': 

^ C = ^ + A- 1 .F C ' (80) 

where the inverse Euler matrix is equal to its transpose A' 1 = A T because of its orthogonality. 
The current sphere's global Z component is given by 

Z c = Zi + (a c + oc,)(sin0ji sin9 c ' coscp c ' + cosOjj cos9 c '). (81) 

Its south pole hits the catch net if Z c = Z net + a c . Solving for the primed azimuthal angle <p c ' at 
which this equation may be satisfied, we have 



Z net + a c ^ _ CQS ^ vC0S ^t 

cos <p c ' = — (82) 

sin0,,sin0 c ' V } 



where the sign of q> c ' is the same as the sign of cp ci \ the angle at the top of the roll corridor. If this 
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equation does not have a real solution in [0 5 7i/2], then the current sphere does not reach the catch 
net in this roll corridor. If (p c ' is outside the range [0, n/2], then units of 2k are either added or 
subtracted from cp c ' until it does lie within that range. Because neither 0jj nor 0 C ' can ever be zero 
or 7C, Eq. 82 is never singular. 

[00267] The preferred implementation also can find which water level the current sphere 
might hit. Water level altitudes between adjacent cylinder walls are generally different as the 
pack is constructed. The only condition under which the current sphere might contact the A>th 
water level Z w (k) is if and only if the distance P c of the current sphere from the global Z axis lies 
between cylinder wall radii Wk-i and Wk at the point of contact. P c is given by 

P c = (X c 2 + Y c 2 ) m (83) 

where X c and Y c can be obtained from Eq. 80 as 

X c = X x : + (a c + a,-)(-cos0y/ cos^ sin# c ' cos#> c ' + sin^y,- sin# c ' sinp c ' 

+ sinOjj cosfyj cos6c) (84a) 

and 

Y c = Yi + (a c + ai)(-cos8ji sin^/ sin^ c ' cos^ c ' - cos^-,- sin^ c ' sin^ c ' 

+ sinOji sirups cos6 c ') (84b) 

[00268] Using the square of P c rather than Eq. 83, the water level checked for contact is 
Z w (k) if W k 2 _ , < PI < W% . Actually, each annular region between the (&-l)-th and k-th cylinder 
walls is checked to see if the roll corridor has a possibility of contacting the water level there. 

[00269] The expression for the <p c ' at which the south pole of the current sphere touches 
the k-th water level is the same as Eq. 82 with Z net replaced by Z w (k). Then the end of the roll 
corridor <p C f is either ±nl2 or a contact point with the catch net or a water level as given by 
Eq. 82. 
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Losing contact with one or both spheres 

[00270] The preferred implementation also can determine the q> c ' at which the current 
sphere will lose contact with either Sphere i or j or both if it encounters nothing else as it 
descends its roll corridor. 

[00271] Of the two spheres, / and j, the one whose center is higher by the subscript u and 
the one whose center is lower is denoted by 7. Their global positions are (X u , Y u , Z u ) and (Xi, Y u 
Z/) and their radii are a u and a/, respectively. We always have that Zu > Z\. Consider the 
unprimed coordinate system whose origin is in the center of the lower sphere and whose axes are 
parallel to the global axes. Let 0 C and (p c be the polar and azimuthal angular position of the 
current sphere in the unprimed local coordinate system centered in the lower sphere. 

[00272] As before, the primed coordinate system is the system obtained by tilting the 
unprimed z axis into the center of the upper sphere and then rotating the xy axes of that tilted 
coordinate system about its z' axis until its x r axis is pointing as steeply upward as possible, i.e., 
has a maximum positive projection onto the unprimed z axis. This rotation is what Eq. 72 does if 
i is replaced by / and j by «. The polar and azimuthal angular positions of the current sphere in 
the tilted coordinate system are 6c' and (pc', respectively where cos0 c ' is given by Eq. 78 with the 
same change in indices mentioned in the previous sentence. 

[00273] The current sphere will lose contact with the upper sphere first. If both lower 
and upper spheres are at the same altitude, it will lose contact with them both simultaneously. 
To find the cpc' position at which the current sphere loses contact with the upper sphere, the 
preferred implementation finds the two longitudinal lines of the lower sphere in the unprimed 
coordinate system along which the current sphere would just nick the upper sphere if it were 
descending along those longitudinal lines. The two (p c at which that happens are the ends of the 
roll corridor at which the current sphere would lose contact with the upper sphere and roll on the 
lower alone. Because the current sphere is rolling down only one side of the roll corridor, only 
that side's (p c is of interest. 

[00274] In the unprimed local coordinate system, the current sphere is at 

/\ /\ 

r c = (a c + a,)(i sin 6 C cos G c + j sin 9 C sin q> c + k cos 9 C ) (85) 



66 



and the upper sphere is at 

r u = r u (i sin 6 U cos 0 U + j sin 9 U sin + k cos 9 U ) . (86) 
[00275] Consider the equation 

\r c -r u \ = a c + a u . (87) 

[00276] If both sides are squared and written explicitly, we obtain 

[(a c + aj) 2 sin0 c cos(p c - Xuf + + a/) 2 sin9 c sincp c - ^ u ] 2 + 

[(a c + a/) 2 cos0 c - z u ] 2 = (a c + a M ) 2 (88a) 

which can be rewritten 

r 2 +{a c + ty) 2 -{a c + a u ) 2 . _ . _ . _ 

— = x u sin # c cos <p c + j/ B sm 0 C sin (p c + z u cos 0 C . (88b) 

2(a c + a } ) 

[00277] The components x u , y U9 and z u are given in Eq. 86. If they are substituted into 
Eq. 88b and both sides are divided by r U9 we obtain 



+(a c + a l ) 2 -(a c + a u ) 2 . _ . _ 

= sin 0 U cos <p u sin 6 C cos q> c + 

2(a c + a,K * " ^ (g9) 

sin 6 U sin sin 0 C sin <p c + cos # w cos # c 



where we recognize the left-hand side as cos0 c ', the polar angle in the primed coordinate system 
which is constant as the current sphere descends the roll corridor. Now Eq. 89 can be written as 

cos9 c ' = sinG u sin0 c (coscp u coscp c + sin<p u sincp c ) + cos0 u cos0 c (90a) 
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which through a trigonometric identity can be written 



cosO c ' = sin0 u sin9 c cos((pc - (p u ) + cos9 u cos0 c . (90b) 

[00278] We wish to solve for 0 C from this equation. For a given cp c , the 0 C solution has 
three possibilities: 

[00279] 1. There are two 0 C solutions if the current sphere hits the upper 
sphere as it rolls down the (p c longitudinal line of the lower sphere - one solution where it 
touches the upper sphere from above and the other where it touches from below; 

[00280] 2. There is one solution (the two solutions having become 
degenerate) where the current sphere just nicks the upper sphere while rolling along the cp c 
longitudinal line; and 

[00281] 3. There are no real solutions, meaning the current sphere 
altogether misses the upper sphere while rolling along the <p c longitudinal line. 

[00282] The second possibility is the one of current interest because that is the point at 
which the current sphere can lose contact with the upper sphere and roll on the longitudinal line 
of the lower sphere. We will call that point the "departure point." Therefore, the preferred 
implementation first finds the angular position (0 C , (p c ) of the departure point and then finds the 
azimuthal angle <p c ' in the primed coordinate system to which this angular position corresponds. 

[00283] To solve for & C9 given a particular <p c , the preferred implementation recasts 
Eq. 90b into the form 



cos 6 ' 

■ = sin a sin 0 C + cos a cos 9 C (91) 



yjsin 2 G u cos 2 (<p c - <p u ) + cos 2 9 U 



where 



sin a = Bmg.co8te-g,,) 

ysin 2 9 U cos 2 {(p c - (p u ) + cos 2 8 U 
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and 



cos a = 



cos 9 U 



yjsin 2 9 U cos 2 {<p c - <p u ) + cos 2 6 U 



(92b) 



[00284] Applying our trigonometric identity to the right-hand side of Eq. 91, we obtain 



cos(0 c -a)- 



cos 6 r ' 



^/sin 2 6 U cos 2 (<p c - #> u ) + cos 2 # u 



(93) 



which gives the solution 



0 c = a + cos" 1 



cos # c ! 

^/sin 2 # u cos 2 - + cos 2 0 U 



(94) 



[00285] The arccosine generally has two distinct solutions: (1) the principal value 
which lies in the interval [0, 7t] and which is denoted by capitalizing the function as Cos" 1 and 
(2) the negative of the principal value. We are interested in the case where these two solutions 
are degenerate. That can only be so if the arccosine and its negative are the same value, i.e., 
zero. Therefore 0 C = a at the departure point, and the argument of the arccosine in Eq. 94 must 
be unity: 



cos e: 



^/sin 2 6 U cos 2 (<p c - <p u ) + cos 2 6 U 



= 1 



(95) 



[00286] This equation is the condition used in the preferred implementation to find the 
unprimed azimuthal position of the departure point. Squaring both sides gives 

cos 2 9 c '= sin 2 9 u cos 2 (cp c - cp u ) + cos 2 9 u (96) 
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which can be q> c to yield 



<Pc = (p u + COS 



/cos 2 6c -cos 1 0 U 
' sin 2 0 u 



(97) 



[00287] The square root has a positive sign because the difference between (p c and <p u 
can be no greater in magnitude than nil. Again, the arccosine has two solutions: the principal 
value Cos" 1 and its negative -Cos" 1 . The current sphere would lose contact with the upper sphere 
if it is at azimuthal position <p u ± Cos" 1 . Its true position depends on where it was initially. If its 
initial position <p C j > q> u , then the principal value Cos -1 is used. If its initial position (p C j < (p u , then 
the negative of the principal value is used -Cos" 1 . In the unlikely case that the initial position 
equals <p u , then the current sphere is balanced on the upper and lower sphere and a direction is 
randomly chosen by the random number generator. The preferred implementation using these 
relationships is to determine where contact is lost. 

[00288] The preferred implementation also determines (0 C , cp c ) of the departure point. It 
then finds the departure point in the primed coordinate system by use of the Euler transformation 
of Eqs. 72 and 76. 

[00289] The first component of the primed vector is given by 

sin9 c ' coscpc' = -cos0 u sin0 c (coscp u coscpc + sincp u sincp c ) + sin0 u cos0 c (98) 

where the expression in parentheses can be written cos(cp c - (p u ). Eq. 97 can substitute for cos(cp c 
- (p u )- Also because 0 C is equal to a of Eqs. 92a and 92b, we have 



. ^ Jsin 2 6>„-sin 2 0 c ' 

sin 0 C = ± (99a) 



cos Oc 



cos 9 U 

cos0 c = (99b) 
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[00290] Substituting Eqs. 97 and 99a/99b into 98, we obtain for the end of the roll 
corridor due to departure 



, cos# w sin# c ' tm0 c ' „^ 

cos «v = —f — fr = — f • ( 10 °) 

sin0 B cos0 c tan# u 

[00291] The preferred implementation uses this relationship to obtain 9 C and cp c . Notice 
that the current sphere cannot be in stable contact with two pack spheres unless 9 c f < 8 U . 
Therefore (pcf always lies between 0 and 7t/2. Also, remember that cosG c ' is given by the 
left-hand side of Eq. 89 and that sin#c' = +[1 - cos 2 6 z '] m . 

[00292] The second component of the primed vector is given by 

sin9 c ' simpc' = (sin(p u coscp c - coscp u sin<p c ) sin9 c (101) 
and the same set of substitutions leads to 



, Jsin 2 0 w -sin 2 0 c ' 
sm <P<f = " t . „ • (1 02) 



[00293] The (p cf ' defined by Eqs. 100 and 102 is the departure point and is the end of the 
roll corridor. At this position, the current sphere loses contact with the upper sphere (or both if 
their centers are the same altitude). These relationships also are used in the preferred 
implementation to determine 9 C and (p c . 

Hitting the cylinder wall 
[00294] The cylinder wall for the current sphere is located at global radius W. The 
preferred implementation can find whether the roll corridor of the current sphere on the two pack 
spheres is such that the current sphere's center can intersect the cylinder wall. If it does, then the 
cp c ' of intersection is found. 
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[00295] We need the global position of the current sphere as a function of cp c '. Again, 
the global position of the current sphere is 

Rc = Ri + r e (103) 

where J?,- is the global position of one of the contact spheres and r c is the current sphere's 
position in the unprimed local coordinate system. (Remember the unprimed coordinate system is 
centered in Sphere i and has axes parallel to the global axes.) The unprimed vector is given in 
terms of the primed vector (that contains cp c ') by the inverse Euler matrix equation 



F c = A -1 • r c f (104) 

where A' 1 is the transpose of the Euler matrix in Eq. 72 and r c ' = (x c \ y c \ z c ') is given by 

*c = (o-c + a t ) sin9 c ' cos(p c ' (105a) 

Vc " («c + «i) sin9 c ' sin cp c ' (105b) 

= (a c + a,-) cos0 c '. (105c) 



[00296] The preferred implementation therefore finds the global position R c = (X C9 Y c , 
Z c ) in terms of the single variable <p c ' describing the descent down the roll corridor. If the lateral 
global position of the current sphere reaches the cylinder wall radius W for a value of <p c ' in the 
interval (<p C j' <p C f ), then the current sphere can reach the cylinder wall if it does not encounter 
another object first. In other words, if there is a solution of 

X C 2 + Y C 2 =W 2 (106) 

for (p c ' in the interval ((p c i' 5 q> C f )> th en the current sphere reaches the wall. Written explicitly, 
Eq. 106 has the form 

A 0 + A j cos(p c ' + A 2 cos 2 (p c ' + B } sin(p c ' + B 2 sin 2 cp c ' = 0 (107) 
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where 



A 0 = (Xj 2 + Yl - W i /(a e + a,) + sin8 y , cos0 c ' [2X ( coscpy,- + 





2Yi sincf*,,- + (a c + a,) sinGy, cos9 c '] 


(108a) 


A, 


= -2 sin9 c ' COS0// [Xi costy, + 7, sincpy,- + (a c + a,) sinG/, cos9 c '] 


(108b) 


A 2 


= (a c + a,) cos Qji sin 0 C ' 


(108c) 


B, 


= 2 sin9 c ' [Xi sincpy, - 7,- cos(p yi ] 


(108d) 


B 2 


= (a c + a,) sin 2 0 c '. 


(108e) 



[00297] Eq. 107 can be converted into a quartic equation in sinq> c ' by changing the 
cos 2 (p c ' term into 1 - sin 2 (p c ' and the cos(p c ' term into +[1 - sin 2 q> c '] 1/2 . (The plus sign is used only 
because <p c r cannot exceed ±tt/2.) We then have 

(A 0 + A 2 ) + B } sin9 c ' + (B 2 - A 2 ) sinV c ' = -Ai[\ - sin 2 (p c '] 1/2 . (109) 

[00298] Squaring both sides gives 

(A 2 - B 2 ) 2 sinV' - 2Bj(A 2 - B 2 ) sin 3 ^' + (Aj 2 - 2AqA 2 - 2A 2 2 + B } 2 + 2^o5 2 + 

2^2) sin 2 ^ c ' + 2(^ 0 + A 2 )Bj sin^ r + (^ 0 - A } + ^ 2 )(^I 0 + ^ 7 + ^ 2 ) = 0 (110) 

[00299] If this equation has any real roots in the interval (<p ci ' 9 (p cf '), then the smallest in 
magnitude of these that satisfies the original Eq. 107 is the point at which the roll corridor 
intersects the cylinder wall. The preferred implementation can check the roots in the original 
equation because squaring can give extraneous roots. 

Hitting one or more spheres simultaneously 
[00300] The preferred implementation has now defined a roll corridor in the primed 
coordinate system beginning at <p c / and ending at (p c f. The end of the roll corridor is where the 
current sphere either lost contact with one or both of its contact spheres, when its south pole 
encountered either a catch net or a water level, or when its center encountered its cylinder wall. 
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The preferred implementation can check all spheres that can intersect this roll corridor to find the 
sphere(s) encountered first as the current sphere descends along its roll corridor, i.e., as <p C ' moves 
from (pd to <p c f. 

[00301] The global position of the current sphere as a function of q> c f is given by Eqs. 80, 
81, and 84. If some other pack spheres, say the &-th sphere is to intersect the roll corridor, then 
the following equation must have a solution in the range {(p ci \ (p c f)\ 

R c {<p c ')-R k \ = a c + a k . (Ill) 
[00302] Squaring this equation gives 

R 2 C -2R C -R k + R 2 k =(a c + a k ) 2 (112) 
which becomes the following equation in <p c ': 

A cos(p c ' + B sm<p c ' = C (113) 

where 

A = -X ik cos tpji cos Oji sin d c ' - Y ik cos sin (p^ sin 0 C ' + Z ik sin 9 C ' sin 0p (114) 
B = X ik sin (pji sin 9 C ' - Y ik cos (pji sin 0c f (11 5) 

C = -X ik cos (pji cos 8 C ' sin dji - Y ik cos d c ' sin ^ sin Q j{ - Z,* cos # c ' cos 0, 7 
t (2a c + a, + g t )(fl ft - a f ) - ^ 

+ . (116) 

2(a c + a,) 

[00303] The global coordinates of Sphere Z: relative to Sphere i are given by 

Xik - X/c, Y ik = Yi - Yk, Zik = Zi - Z k , (1 17a) 



and the square of the distance between the two sphere centers is 
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Rik - Xik + Yfr + 



(117b) 



[00304] The solution procedure of Eq. 113 is similar to that given in Section VII.D.7. 
The preferred implementation finds <p c ' to be 



where 



C 

(p c ' = a± cos" 1 , (118a) 

<JA 2 + B 2 



a = tan" 1 -. (118b) 
A 



[00305] If any of the <p c f solutions are outside the range [0, In], units of 2k are added or 
subtracted until they are inside that range. Then if any solutions lie within {(p ci \ <p c f) the smallest 
in magnitude of these is the position of first contact between the current sphere and Sphere k. If 
there is no solution, then Sphere k intersects no part of the roll corridor. 

[00306] The preferred implementation has now performed an exact solution to find 
whether and where Sphere k intersects the current sphere's roll corridor. In most packs, the vast 
majority of the pack spheres are out of range of the roll corridor and so do not need the detailed 
calculations of Eqs. 1 18a and 1 18b. There are a few much simpler checks that can be performed 
by the preferred implementation to immediately determine whether a pack sphere is out of range 
and so is not a contact candidate. 

[00307] First, the roll corridor's highest global altitude is at <p ci ' which, from Eq. 81 is 
found to be 

Zcjnax = Zi + (a c + a,-)(sin0ji sin0 c ' cos(p ci ' + cosGji cos0 c '). (1 19) 

Therefore, the k-th pack sphere whose center lies a distance a c + ak above Z Cimax cannot overlap 

the roll corridor. Likewise, the lowest global altitude is at cp c f and is given by 
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Zcjnin = Zi + («c + a/)(sin0ji sin9 c ' cos(p c f + cosGji cos0 c '). (120) 

[00308] The A>th pack sphere whose center lies a distance a c + ak below Z C)min also 
cannot overlap the roll corridor. Therefore, the preferred implementation can immediately 
exclude from consideration all spheres whose global altitudes Z* satisfy the following: 

A>th sphere does not overlap if 

Z k > Z c , max + a c + a k (121a) 



or 



Zk<Z cMn -a c -a k . (121b) 

[00309] Roll corridor extrema in the X and Y directions do not necessarily lie at the 
beginning and end of the roll corridor. Eqs. 84a and 84b give the global Xand Y coordinates of 
the current sphere as a function of (p c \ If we differentiate X c {(p c ') with respect to (p c * and set the 
derivative equal to zero, we have an equation whose solution gives us the extrema in X, If either 
of those extrema lie in {(p C i\ <p c f) then that position is either the minimum or maximum X value 
the current sphere assumes as it descends the roll corridor. The same is true for the derivative of 
Y((pc) with respect to <p c '. The derivatives followed by the extrema solutions are 



dX 

c - = (a c + a,Xcos Oji cos q> fl sin 0 C 1 sin <p c ' + sin (p }i sin 6 C 1 cos <p c ') = 0 (1 22) 



d<p t 



t»fW~*Sf* (123) 

cos 8 ji 



for the two A" extrema (the above equation has two solutions for <p cx '), and 
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= (a c + a,)(cos 6ji sin (p yt sin 6 C ' sin cp c ' - cos <p ;i sin 9 C ' cos <p c ') = 0 (1 24) 



I = _cot^ (125) 



COS 



for the two Y extrema. 

[00310] Finally then, the X and Y extrema of the roll corridor in the global coordinates 

are 

X c<min = MW[Xc{(p ci '), X c ((p c /), X c (<p cx ')] (1 26a) 

X^at = MAX[* c (p c /), X c (<p cf % X c (<p cx ')] (126b) 

with the X c ((p cx ') omitted if neither yd lies in (<p ci ', <p c f). Similarly, 

Y c , min = MM[Y c (<p ci '), Y c {(p c f), 7 C «)] (127a) 
Y Ctmax = MAX[Y c (<p ci % Y c ((p c /), Y c (<p c ;)], (127b) 

again with the Y c {(p C y) omitted if neither <p cy ' lies in (<p C i', (p c f)- 

[00311] Then the &-th sphere cannot overlap the roll corridor if 

X k <X cMn -a c -a k (128a) 



or 



X k >X Cimax + a c + a k (128b) 



or 



Y k < Y c>min - a c a k (128c) 
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or 



Y k >Y Citnax + a c +a k . (128d) 

[00312] Another discriminator that can be used by the preferred implementation to 
immediately exclude a large number of pack spheres from detailed examination is to consider 
whether the roll corridor is far enough from the central Z axis that it cannot overlap all the 
spheres within a given inner cylinder. To find this, the preferred implementation uses the closest 
radial approach P c ,mm of the roll corridor to the central axis. The square of the radial distance 
[Pc(<Pc)] 2 is given by 

[PciVc')? = [XM? + [Y c (<p c ')] 2 (129) 
where we obtain the explicit expression by substituting Eqs. 84a and 84b into Eq. 129: 

[Pc(<Pc)] = \Xi + («c + ai)(-cos#y,- cos^/7 sin# c ' cos^ c ' + sin^? sin# c ' sin#> c ' 

+ sin0y/ coscpji cos 0 c r )] 2 
+ [Yj + (a c + ai)(-cos6ji sin^ sin# c ' cos^ c ' - cos^y,- sin# c ' sin^ c ' 

+ sinfy sin^y,- cos<9 c ')] 2 . (130) 

[00313] Setting the derivative of [P c (<Pc)] 2 with respect to <p c ' equal to zero yields the 
equation whose solution <p cp ' gives the closest approach of the roll corridor to the global Z axis: 

[Xi + (a c + ai)(-cos8ji costpji sin#c' cos^ c r + sirupji sin# c ' sm<p c ' + 

sin^v cos<pji cos# c ')K cos #// cos^ y7 sin# c ' sm<p c ' + sin^v sin# c ' cos<p c ') + 

[7/ + (a c + ^(-cos^y, sinfyi sm6 c ' cos<p c ' - coscpp sin#c' s\n(p c ' + 

smOji sm(pji cos0 c ')]( cos ^7 sin^y, sin#c' sin^ c ' - cos^y,- sin# c ' cos^ c ') = 0.(131) 

[00314] This equation is of the form 
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Aq +Ai cos(pc + A 2 cos 2 (pc + Bj sm<p c r + B2 sin 2 ^ c ' = 0 



(132) 



where in this case, the coefficients are given by 

Aq ^Xi + Y/ 2 + 2(a c + a t ) cos0 c '(A!/ sinfy- cos^-,- + Y t sin^- sin^-,-) + 



(a c + a/) 2 cos 2 0 c ' sin 2 0 y7 (133a) 

A j = -2(a c + aft cosOji sin# c ' [Xi cos^- + 7/ sin^-,- + 

(a c + A/) cos# c ' sin^/] (133b) 

^2 = (fl c + 0/) 2 cos 2 ^/ sin 2 ^' (133c) 

5y = 2(a c + a,-) sin# c ' {sin0y ( - [Xi + (a c + a,-) cos# c ' sinfy cos^] 

- cosd c '[Yi + (a c + a/) cos# c ' sin^/ sin^,-]} (133d) 

#2 = (a c + af sin 2 0 c ' (133e) 



[00315] Eq. 132 can be converted into a quartic equation by changing the cosines into 
sines in the same fashion as led up to Eq. 110: 

(A 2 - B 2 ) 2 sinV - 25/0*2 - B 2 ) sin V + (A, 2 - 2AoA 2 - 2A 2 2 + 
Bj 2 + 2AoB 2 + 2A 2 B 2 ) sinV' + 2{A 0 + A 2 )Bj sin<p c ' + 
(A 0 -A; + A 2 )(A 0 + A } + = 0. (1 34) 

[00316] There are four sin <p c ' solutions to this quartic equation and for a given sin <p c ' 
there are two <p c ' solutions. Complex solutions and real solutions lying outside [1, 1] are 
immediately discarded in the preferred implementation. Furthermore, we are only interested in 
those solutions with the same sign as sinp c / and that lie in {(p ci \ (p c f) because that is the direction 
the current sphere will roll. The preferred implementation checks any remaining solution <p cp 9 by 
ensuring it satisfies the original Eq. 132. If there is no solution in (<p ci \ q> c f) then one or the other 
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of the end points is the closest approach to the global Z axis and both must be substituted into 
Eq. 130 to find the smaller P c . 

[00317] Once the P c , m in of closest approach is found, the preferred implementation 
checks to see if there are any mode's cylinder radii Wk that lie completely inside P c , m in by more 
than the sum of the mode sphere's radius plus the current sphere's radius. If 

Pc>W k + aj+a c (135) 

then none of the Mode k particles can touch the current sphere as it descends its current roll 
corridor and so may immediately be excluded as candidates for contact. 

[00318] All spheres whose north poles have sunk beneath the catch net or water level are 
not candidates for contact with the current sphere. The preferred implementation continually or 
periodically monitors the lowest index number of the sphere whose north pole is still above the 
catch net so that no time is wasted looking at spheres with lower index numbers as candidates for 
contact. 

Touching one sphere and the cylinder wall 

[00319] Remember that in the preferred implementation it is the center of the current 
sphere constrained by the cylinder wall, not its edge. When in stable contact with both a pack 
sphere and the cylinder wall, the roll corridor lies on the wall and descends toward the equatorial 
plane of the contacted pack sphere. If the contacted sphere (denote it Sphere i) is as large as the 
cylinder, then the current sphere moves toward the lowest point on a roll corridor that spans the 
entire 2n radians of the cylinder. 

[00320] The roll corridor will be denoted by the global azimuthal position & c . Its initial 
position is given by 

0 cO = tm\Y c /X c ). (136) 

[00321] The preferred implementation uses the range [0, 2n] rather than [-7C, n] for the 
range of the global azimuthal position. The global azimuthal position of the Sphere i is 
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&i = tan 1 {YilXii 



(137) 



which also lies in [0, In]. If & c o > &i the current sphere rolls toward higher 0 If 0 cO < <P„ it 
rolls toward lower <2>. In the unlikely event that & c0 = 0i, the random number generator chooses 
a roll direction. 

[00322] If it encounters nothing else, the current sphere will roll until it reaches: 

[00323] 1. the equatorial plane of Sphere i where it will fall away from 

Sphere z; or 

[00324] 2. the lowest azimuthal point in its cylinder permitted by Sphere i. 
[00325] The azimuthal position at which it will reach the equatorial plane of Sphere / is 
found by the law of cosines as 

0 cf = 0 i ±A0 c (138) 



where 



A<P c = cos _1 — v J 



2W C P { 



(139) 



[00326] In this expression, W c is the radius of the cylinder wall, P ( = ^A? + Yf is the 
radial position of Sphere / from the global Z axis, and a c and a t are the radii of the current sphere 
and Sphere i 9 respectively. The plus sign is used in Eq. 138 if the current sphere is rolling toward 
increasing 0, the minus sign is used if the current sphere is rolling toward decreasing 0. If the 
absolute value of the quantity in the parentheses of Eq. 139 is greater than unity, then there is no 
solution and the current sphere simply rolls around the perimeter of the cylinder until it reaches 
the lowest point it can. This lowest point is opposite the azimuthal position of Sphere /. Hence 
if 
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>1 



(140) 



2W C P { 



then 



A0 C = ±7T . 



(141) 



[00327] In the unlikely event that P t = 0, i.e., Sphere i is centered on the global Z axis, 
then there is no downward direction. The roll corridor is horizontal and the current sphere is 
placed where it is. 

[00328] In this analysis, if the addition or subtraction of A0 C causes 0 c /to lie outside [0, 
2n] y then units of 2% are added or subtracted until & c f again lies in [0, In}. 

[00329] The preferred implementation also determines what placed pack spheres might 
intersect the roll corridor in the range (0 cO , 0 C /). 

[00330] The current sphere will impact a placed pack sphere, say the >th sphere, if there 
is a solution to 



anywhere in the range (& c0 , <P c f). The Cartesian components of the vector R c (O c ) = X c , Y c , Z c are 
given by 



[00331] The Z c component can be found by inserting the X c and Y c component into the 
condition 



\R c (O c )-Rj =a c + aj 



(142) 



X C (0 C ) = W C COS& C 



(143) 



Y c (0 c ) = W c sm0 c 



(144) 
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\R C -Ri =a c + ai . (145) 

[00332] Squaring both sides, we obtain 

W c 2 + 2c + i?, 2 - 2W c (Xi cos0 c + Yi sin0 c ) - 2Z,Z C = (a c + a,) 2 (146) 

which is quadratic in Z c . One of the solutions is when the current sphere is above Sphere i and 
the other when it is below. In the preferred implementation, we are only interested in the 
solution above Sphere i, which is 

Z C (0 C ) = Zi + [Zf + (a c + a,) 2 - W c 2 - R 2 + 2W c (X t cos0 c + Yi sin<Z> c )] 1/2 . (147) 

[00333] Now, inserting Eqs. 143, 144, 147 in Eq. 142, we obtain an equation to solve for 
the & c at which the current sphere contacts the j-th pack sphere: 

A 0 + A, cos& c + B, sin& c + (A 2 cos& c + B 2 sin<£ c ) 2 = 0 (148) 

where the coefficients are given by 

A 0 = G 2 - 4(Zi - Zj) 2 [Z? - W c 2 - R? + (a c + a,) 2 ] (149a) 

A, = AGWdXi - Xj) - 8(Z, - Zj) 2 W c Xi (149b) 
A 2 = 2W c (X i -X j ) (149c) 

B, = AGW c (Yi - Yj) - 8(Z, - Zjfwji (149d) 
B 2 = 2W c (Y i -Y J ) (149e) 

and G is given by 
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G = (a,- - aj)(2a c + a, + aj) - R? + R? + 2Zf - IZiZj . 



(150) 



[00334] Eq. 148 can be converted into a quartic polynomial in cos0 c : 



{A 2 + B 2 2 ) 2 cos 4 0 c + 2{A } A 2 + 2A 2 B } B 2 - AjB 2 2 ) cos 3 0 c + 
(A j 2 + 2AqA 2 2 + Bj 2 - Z4oB 2 2 - 2A 2 2 B 2 2 - 25/) cos 2 tf> c + 
2(AoAj - 2A 2 BiB 2 + A } B 2 ) cos& c + 

(^-5y+£/)(^ + 5y+3/) = 0 



(151) 



whose solutions are candidates for contact of the y-th pack sphere. The preferred implementation 
uses this expression as part of its determination. There are four solutions of the quartic equation 
and each cos& c solution has two 0 C solutions — eight in all. Of the four quartic solutions, the 
preferred implementation can eliminate complex solutions or solutions whose magnitude is 
greater than unity. Of the remaining solutions, if there are any that satisfy the original Eq. 148 
and that lie in the range of the roll corridor (0 cO , 0 cf ) y they are contact points. Of the contact 
point(s), the value of 0 C that lies closest to 0 cO is the point at which the current sphere touches 
the 7-th sphere. We will denote that point by & cj . If there are no solutions that satisfy the above 
conditions, the current sphere cannot contact the 7-th sphere. 

[00335] As the preferred implementation checks the candidate spheres in the pack for 
contact with the current sphere, that sphere with 0 cj closest to 0 cO , if any, is the next sphere 
contacted. In the unlikely case that two or more pack spheres have the same 0 CJ , then it checks 
the stability of each trio of contacts to see if the current sphere is to be placed there. If it is not 
stable, then the stability routine determines which sphere(s) the current sphere will roll away 
from. 

[00336] As an efficiency measure, the preferred implementation can eliminate from 
consideration any modes or categories that are sufficiently far inside the cylinder wall of the 
current sphere that they cannot touch it. All mode k spheres are out of range if 



W k + a k < W c -a c . 



(152) 
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[00337] Any other sphere j is a candidate for being in range only if 



W c -a c -aj<Pj< W c + a c + aj 



(153) 



where Pj is the global radial distance of Sphere j from the Z axis. 

THE SUBROUTINE UPDATE 

[00338] The UPDATE subroutine of the preferred implementation updates all 
parameters and arrays that change as a result of the placement of the current particle. These 
changes include the following. 

1. The Pool Switch 

[00339] If the pool switch is true, then the placement of the current particle means the 
bin from which it was chosen should be diminished by one. This adjustment is performed in 
Subroutine UPDATE. 

2. The Catch Net 

[00340] The catch net belonging to the mode, submode or category current sphere also 
generally ratchets upward. The prescription for determining how much it ascends will now be 
addressed. 

[00341] A sense of pack surface must be established. With a variety of particles being 
dropped into a variety of cylinders, the pack surface must be defined. This can be done in a 
number of ways, as noted above. For purposes of the preferred implementation, let us consider 
the altitude of the south poles of the last m placed spheres of any and all modes. As explained 
above, the average of those altitudes is defined for present purposes as the pack "surface": 



'surface 



1 m 

= -H(Zj-aj) 
m 7-i 



(154) 
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where m is the number of most recently placed pack sphere equatorial cross-sections that would 
just cover the outermost cylinder n with cross-sectional area nW n 2 . The value for m is found 
from the inequality 



m-\ m 



Sa/ <w 2 n <Za/ . 



(155) 



7=1 >1 



[00342] If the size of the last placed set of spheres is very large so that m is very small, 
then we have a minimum required value for m which can be user-defined but for which the 
default value in this implementation is m min = 12. Only when the pack is just starting so that 12 
particles have not yet been dropped is m smaller for this implementation. Preplaced particles do 
not count in assessing the pack "surface." 

[00343] The k-th mode catch net is then defined as 



where y is the number of A:-th mode radii beneath the surface that the A>th mode catch net lies. 
Here it is a user-defined number which has a default of y = 2. 

[00344] Subroutine UPDATE readjusts all the above parameters so that the catch net for 
every mode has a proper value after the current sphere is placed. 

[00345] When the pack is still small enough that an insufficient number of particles have 
been placed to satisfactorily define a surface, the catch net may be ratcheted up a little after the 
placement of each particle so that the particles on the bottom of the cylinders do not all have the 
same altitude. If the switch ROUGHNET is true, then the current sphere's catch net moves 
upward by a small random amount 



where r is a uniformly generated random number in [0, 1] and N cs is the total number of sphere 
cross sections required to cover all the cylinders: 



Znet(Jt) 2 surface " J&k 



(156) 



ra min IN c 



cs 



(157) 
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(158) 



where n is the number of modes, Wk is the cylinder wall of the A>th mode, and Wo = 0. If 
ROUGHNET is false, then the catch net stays where it is until enough particles have been 
defined to make a surface and until the surface is more than Z sur f ace - yak above the starting place. 

[00346] The starting place for the k-th mode catch net is Z init - a*. A £-th sphere radius is 
subtracted from the overall pack starting altitude Z/ m , because the catch net stops south poles and 
Z init is the starting place for sphere centers. The default for Z init is zero in this implementation. 

The Water Levels 

[00347] The water level also is modified when the current sphere is placed according to 
the preferred implementation. Recall that the water level only exists beyond the innermost 
cylinder. Its altitude tracks the pack surface and is at a user-specified distance p in terms of the 
smallest sphere radius a min beneath the surface defined by Eqs. 154 and 155. The default for p is 
unity. Hence 

Z w — Z sur f ace - p<2 /m/I . (159) 

If P is negative, the water level lies above the pack surface. 

[00348] The smallest cylinder in this implementation has no water level, only a catch 
net. It does not need a water level because water levels are designed to simulate the surface of 
the innermost cylinder and the innermost cylinder has its own particle surface potentially 
containing all the particles of the pack. 

Current Sphere History Parameters 

[00349] The preferred implementation uses are parameters that monitor the number of 
times the current sphere tunneled, the number of distinct roll corridors it underwent (not counting 
free-fall drops), and the number of free-fall drops it underwent. All these are reset to zero in this 
implementation. 
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[00350] The arrays indicating which pack spheres the current sphere is touching and 
which spheres it just rolled away from are also zeroed out along with the total number of such 
spheres. 

[00351] The array keeping track of the lowest index number above catch net for each 
mode is adjusted if necessary. 

[00352] The array monitoring the index numbers of all pack spheres with exposed north 
poles also is updated if the current sphere covered one or more exposed north poles. 

[00353] In addition, the pack sphere position arrays have the current sphere's (X, Y, Z) 
position added to them along with the array specifying its mode number. The counter 
monitoring the total number of spheres placed so far in the pack is also incremented by one. 

THE SUBROUTINE OUTSTAT 

[00354] There are a number of important pack statistics and features that the user may 
need or want to calculate. In this implementation, subroutine OUTSTAT calculates the output 
statistics that the user specifies. These output statistics include: 
[00355] A. The volume fraction of the pack. 

[00356] B. The radial distribution functions gij(r)dr which is the mean 
number of Mode j particle centers that lie between r and r + dr from a Mode / particle center. 

[00357] C. Coordination numbers Qy which represent the mean number of 
Mode j particles touching a Mode / particle. 

[00358] D. Detailed pack statistics for each particle mode such as that 
mode's minimum and maximum positions, number of its particles in the pack, radial and axial 
density functions for each requested mode, as well as other information. Each of these output 
data will now be addressed. 

Finding the Volume Fraction 

[00359] To determine the packing fraction, the program finds the total volume of the 
innermost cylinder occupied by spheres. Recall that the cylinder walls are boundaries to the 
centers of the spheres, not their edges. Hence, many spheres will be overlapping the cylinder 
walls. The preferred implementation has a general expression for determining how much 
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volume of a sphere of radius a s located at position (x s , y S9 z s ) lies inside a cylinder which is 
concentric to the z axis with radius a c , bottom cylinder altitude zt and top cylinder altitude z t . 
The sphere volume inside the cylinder will be denoted V in . Here the sphere and cylinder 
geometric parameters have no restrictions except: 

a s >0 9 a c >0 9 and z t > z^ (160) 

[00360] The radial distance of the sphere center from the z axis is defined here by 

p 5 = [^ 2 +^ 2 ] 1/2 . (161) 

[00361] There are four possible cases for this overlap problem. In Case 1, the sphere is 
completely outside the cylinder so V in = 0. In Case 2, the sphere is completely inside the 
cylinder so V in = (4/3)7ta s 3 . In Case 3, the cylinder is completely inside the sphere so V in = 7ta c 2 (z t 
- Zb). In Case 4, the sphere and cylinder intersect each other and the expression for V in is more 
complex. In the preferred implementation, processing is performed to identify and address these 
cases using the following approach. 

Case 1 : Sphere Completely Outside Cylinder 
[00362] The sphere is completely outside the cylinder 

V in = 0 (162) 

only if one of the following conditions are met: 

If the sphere's north pole is lower than the cylinder bottom . . . 

z s + a s <Zb (163a) 

Or if the sphere's south pole is above the cylinder top . . . 

z s -a s >z t (163b) 
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Or if the sphere is more than its radius away from the cylinder side . . . 
p s -as>a c (163c) 



Or if the sphere center is above the cylinder top and is more than its radius away 

from the upper cylinder edge. . . 
(p s - a c f + (z s - z t f > a s 2 and z s > z t (163d) 

Or if the sphere center is below the cylinder bottom and is more than its radius 

away from the lower cylinder edge . . . 
(p s - a c f + (z b - z s f > a s 2 and z s <z b . (1 63e) 

Case 2: Sphere Completely Inside Cylinder 
[00363] The sphere is completely inside the cylinder 

Fa, = (4/3)W (164) 

only if all three of the following conditions are simultaneously met: 

If the sphere's north pole is no higher than the cylinder top . . . 

z s + a s <z t (165a) 

And if the sphere's south pole is no lower than the cylinder bottom . . . 
z s -a s >z b (165b) 

And if the sphere's equator nowhere exceeds the cylinder radius . . . 

p 5 + a s <a c . (165c) 

Case 3: Cylinder Completely Inside Sphere 
[00364] The cylinder is completely inside the sphere 
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V in = ita c \z t - z b ) 



(166) 



if both the following two conditions are met: 



If the upper edge of the cylinder is inside the sphere . . . 

(p s + a c f + (z t - z s f < a s 2 



(167a) 



And if the lower edge of the cylinder is inside the sphere . . . 
(ps + a c ? + (z s - z b ) 2 < a s 2 . 



(167b) 



Case 4: Cylinder and Sphere Intersect 

[00365] There are several subsets of Case 4. To efficiently organize them, we define 
two circles: the sphere circle and the cylinder circle. The sphere circle is the perimeter of the 
sphere's cross section at a given z altitude. Its radius is 



and its radial position p s is given by Eq. 161. Likewise the cylinder circle is the perimeter of the 
cylinder's cross section which has constant radius a c at any z altitude between z b and z t and is 
centered on the z axis. 

[00366] To describe the intersection of sphere and cylinder, four subcases for the 
relationship of the sphere circle and the cylinder circle at a particular altitude z are defined: 

[00367] Subcase 4.1: The two circles do not overlap or one or both are 
nonexistent. 

[00368] Subcase 4.2: The sphere circle is completely enclosed in the 
cylinder circle. 

[00369] Subcase 4.3: The cylinder circle is completely enclosed in the 

sphere circle. 
[00370] Subcase 4.4: The two circles intersect. 



2 2 1/2 

r s = [a s - (z - z s ) ] , z s -a s <z<z s + a s 



(168) 
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[00371] To find V in , the preferred implementation integrates in z the overlap areas of the 
two circles A in from the altitude where one subcase holds to the altitude where another subcase 
holds. 

For Subcase 4.1, A in = 0. (169a) 
For Subcase 4.2, A in = ita s 2 . (169b) 
For Subcase 4.3, A in = na c 2 . (169c) 
For Subcase 4.4, A/ n = a c (a - cosa since) + r s {fi - cos/? shy?) (169d) 
where r s is given by Eq. 168 and 



cosa= P ° +a * r ° (170a) 
2pA 



J[r? - (a c - p 5 f] [(a c + p s f - rj ] 

sina= ^ — v H ' — (170b) 

2p s a c 



cos/3= p ° +r ' a * (171a) 
2p s r s 



• 0 J[^ 2 - (a c - p s f] [(a c + p s f - r? ] 

smfi= — v H ' J LV ysJ — . (171b) 

2p s r s 



[00372] The implementation now constructs a series of altitudes at which the cylinder 
circle and the sphere circle change from one subcase to another in Case 4. 

[00373] A sphere and a cylinder side can intersect each other at either 0, 2, or 4 altitudes. 
Pairs of altitudes may be degenerate (two solutions having the same altitude) if the sphere just 
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nicks the cylinder wall. The preferred implementation determines these altitudes using the 
following equation: 



[r s (z)] 2 = [a c ±p s f (172a) 



or 



a s 2 -(z-z s ) 2 = [a c ±p s ] 2 (172b) 



which, using Eq. 176, has the four solutions 



z = z s ±^af -(a c ±p s f . (173) 



[00374] The + sign in the first ± of Eq. 173 denotes the upper or northern hemisphere 
intersection altitude; the - sign denotes the lower or southern hemisphere intersection altitude. 
The + sign in the second ± of Eq. 173 denotes intersection with the far cylinder wall, on the 
opposite side of the z axis from the sphere's side; the - sign denotes intersection with the near 
cylinder wall, on the same side of the z axis as the sphere. 

[00375] Armed with these four solutions, the implementation can effectively address 
each of the four subcases of Case 4. In this illustrative implementation, it separately addresses 
whether the sphere circle lies in the southern hemisphere of the sphere or whether it is the 
equatorial plane or in the northern hemisphere. The four subcases are considered separately for 
the two hemispheres making eight situations in all. The hemispheres are considered separately 
because if the sphere circle is in the southern hemisphere, it is increasing in size as z increases. 
Hence the altitude at which we will change subcases can be determined. If the sphere circle is 
the equatorial plane or in the northern hemisphere, it is decreasing in size as z increases and 
again the altitude at which subcases will change can be determined. 
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Finding Upper and Lower Boundaries of Subcases for Case 4 

[00376] The logic used in the preferred implementation to find the altitudes at which 
subcases change will now be discribed. Let z/ ovv be the lower altitude of the subcase currently 
being examined and z ni gk be its upper altitude. The function r s (z) is given in Eq. 168. 

[00377] Initially, let z\ ow = zt» the bottom of the cylinder. After each DO WHILE cycle, 
zi ow is set to the previous cycle's z^gh and a new z^gh is sought. To illustrate using pseudocode 
expression: 

V in — 0 ! Initialize V in to zero. 

z/ ovv = z 0 ! Initialize zi ow to the bottom of the cylinder. 

keepGoing = TRUE ! Do Loop flag. 

DO WHILE (keepGoing) 

IF (z s > Ziow) THEN ! Southern hemisphere subcases 

IF (p s - r s (zi ow ) > a c ) THEN ! Subcase 4.1 

^igh = z s - yja? ~(a c - Psf 

IF {zhigh > zd keepGoing = FALSE 

(There is no addition to V in in this subcase.) 
ELSE IF (p s + r s (z low ) < a c ) THEN ! Subcase 4.2 

z high = MTN[z s - y/aj -(a c -p s ) 2 , z„ z s + a $ ] 

Vin = Vin + K{Zhigh ' Zi ow )[a s - Z s + Z s (Zkigh + Z/ow) - 

2 2 
(Z high + ZhighZfow + Z / ow )/3] 

IF (thigh = z, OR z Mgh =z s + a s ) keepGoing = FALSE 
ELSE IF (r s (zi ow ) - p s > a c ) THEN ! Subcase 4.3 

Zhigh = MJN[z s + Jaf - (a c + p s f , z t ] 

Vin Vi n HQ-c (.Zhigh ~ Ziow) 

IF (zhigh = z t ) keepGoing = FALSE 
ELSE IF {{r s {z hw ) > p s AND r s (z, ow ) -p s <a c AND a c < r s (zi ow ) + p s ) 
OR (p s > r s (z, ow ) AND p s - r s (z hw ) < a c AND a c < p s + r s (z low ))) 

THEN ! Subcase 4.4 

IF {a s > a c + p s ) THEN 
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z high = MIN[z, - J a} - {a c + p s f , z t ] 

z M s h 

V in =V in + ]dzA in = V in + W, (174) 

Z low 

where AVj is given in the next section below. 
IF (zhigh = z d keepGoing = FALSE 

ELSE 

z high = MIN[z, + ^af - (a c - p s ) 2 , z t ] 

z hi g h 

V in =V in + jdzA in =V in + AV 2 (175) 

Z low 

where AF2 is given in the next section below. 
END IF ! End two Subcase 4.4 situations 

END IF ! End Subcase 4.4 

ELSE ! Equatorial plane & northern hemisphere subcases 

IF (p s - r s (z low ) > a c ) THEN ! Subcase 4. 1 

keepGoing = FALSE 

(Because the sphere circle is shrinking away from the cylinder 
circle, there is no chance for additional overlap.) 

ELSE TF(p s + r s (z hw ) < a c ) THEN ! Subcase 4.2 

z high = MIN[z,, z s + a s ] 

2 2 

Vin = Vi n + Tl(Zhigh ~ Zi ow ) [a s - Z s + Z s (Zhigh + Zlow) - 

2 2 
{? high + ZhighZfow + Z /ow)/3] 

keepGoing = FALSE 
ELSE IF (r s (z hw ) - p s > a c ) THEN ! Subcase 4.3 



z high = MN[z s + y]af - (a c + p s f , z t ] 

2 

Vin = V in + Jtflc (Zhigh - Zi ow ) 

IF (zhigh = z t ) keepGoing = FALSE 
ELSE IF ((r s (z, ow ) > p s AND r s {z hw ) - p s <a c AND a c < r s {z hw ) + p s ) 
OR (p s > r s (z tow ) AND p s - r s (z hw ) < a c AND a c < p s + r s (z low ))) 

THEN ! Subcase 4.4 
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z high = MJS[z s + yjaj -(a c - p s f , z t ] 
IF(a s >o c + p,)THEN 



high 



V in =V in + jdzA in = V in + AV 3 

Z low 

where AV3 is given in the next section below. 
ELSE 

High 

V in =V in + \dzA in =V in + AV 4 

2 lov 

where AV4 is given in the next section below. 
END IF 



(176) 



(177) 



IF (zhigh - z t ) keepGoing = FALSE 



END IF 



END IF 



END IF 



! End Subcase 4.4 

! End Subcases 

! End hemisphere IF's 



Zlow Zfjigh 

END DO 



! Prepare for next z interval 
! End DO WHILE loop 



Z high 



Evaluation of the Integral j dzA in 
[00378] We now address the integral 



z hi g h 



jdzA;, 



(178) 



where A in is given by Eqs. 169d, 170a, 170b, 171a, and 171b as 
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with 



C = z-z,. (180) 

[00379] This integral will involve elliptic integrals of the first, second, and third kinds 
which for present purposes are defined as in Abramowitz and Stegun. The incomplete elliptic 
integral of the first kind is given by 

F{<p,k)=\de 1 ■ (isi) 

[00380] The complete elliptic integral of the first kind is denoted by K(k) and is equal to 
F{nl2Ji). Similarly, the incomplete elliptic integral of the second kind is given by 

<p 

E{<p, k) = jd0 7l-* 2 sin 2 0 . (182) 



with the complete elliptic integral of the second kind denoted by E(k) and equal to E(nl2Jc). 
Finally, the incomplete elliptic integral of the third kind is given by 



1 

Tl(<p,n,k)= \d0 . . (183) 

0 J {\-nsm 2 e)yj\-k 2 sm 2 0 
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with the complete elliptic integral of the third kind denoted by Tl(n,k) and equal to Yl(it/2,n,k). 
Various efficient routines for the evaluation of these special functions can be found in the 
literature. 

[00381] Focusing on the integral of each of the three terms in Eq. 179, it is possible to 
perform an integration by parts of the first term of Eq. 179 to dispatch the inverse cosine term. 



a 2 \d£ 



cos 



p] + a 2 -a} +C 2 



2p s a c 



= o c 2 ^cos 1 



p) + a 2 -a 2 +£ 2 



+2a2 < 



(184) 



where 



R U2 = J[a? - (a e - p s ) 2 - (*. + p s f - a 2 ] 



(185) 



The right-most integral of Eq. 184 is an elliptic integral. 

[00382] An integration by parts also may be used for beginning the integration of the 
second term of Eq. 179. 



cos 



- 1 



f p\ -a 2 +aj-C 
2p s y/a 2 -C 2 



(2 



C 2 (3a 2 -C 2 )(a 2 + a 2 - p] - C ) 



2\Dl/2 



a} - — Icos 1 



3(a 2 -C 2 )R 
pi -a 2 +a?-C 2 



+ lal{a 2 c -p 2 s )\dS-± 1 



IPsJaj-C 

1(0? ,30/ + I J^-I fl . W -^__J_ 



(186) 
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where the integrals after the second equal sign in Eq. 186 come from partial fraction separation 
of the integrand after the first equal sign. Each of the integrals involves one or more elliptic 
integrals. 

[00383] The third term of Eq. 179 is already in a form amenable to solution by elliptic 
integrals: 



4k* i/: 



(187) 



[00384] Combining these three expressions in Eqs. 184, 186, and 187 yields the general 
expression for the Abused in Eqs. 174 to 177: 



AV = a 2 ^ 2 )cos' 1 
f 



( „2 



pt + a 2 - a 2 + C 2 



cos 



V 



p j - a 2 + a 2 - C 
2p s ^a 2 ~C 



2psa c 

?2 







( 1 ^ 






J 




K 3 J 



+ -af(a 2 -p 2 )L 0 + 



f5 -al -a] + \p 2 )l 2 + Il 4 - \a' s {a 2 c -p 2 )P-±-E 



(188) 



where 



(189) 



(190) 



(191) 
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and 



E= \d£R vl (193) 

with R m given by Eq. 185. The explicit solutions to the L, P, and E integrals of Eqs. 189 to 193 
depend on the values of C;> 0, a s , and (a c - p s ). These expressions are used in the preferred 
implementation. 

Explicit solution for interior volume F/„ when a s > a c - p s 
[00385] In the case where a s > (a c + p s ) , the preferred implementation follows the C 
integration up the z axis, the southern hemisphere's sphere circle eventually circumscribes the 
cylinder circle unless it hits the top of the cylinder first. The limits of integration are therefore 
given by 

(j = MAX[B,-(z t -z s )] and 6 = MIN[^ 5 -(z b - z s )} (194) 

where A and B are given below. If a 2 > (a c + p s ) 2 , then it follows that a 2 > (a c - p s ) 2 . Therefore, 
we define 



A = ^a] -{a c -p s f 



and 



B = ^a 2 s -(a c + p s ) 2 



(195a) 



(195b) 
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and note that A > B > 0. Then 



R^ = ^-C){C~B 2 ). (196) 



[00386] The integrals for this form of R in are readily found in references such as 
Gradshteyn and Ryzhik. Define 



A ' =sin " 1 ,RH^T > (197a) 



A 2 -B 2 ' 



A2 = sin vS^' (197b) 



and 



1=*— A • (197c) 

A 



[00387] Then, recalling the definitions for the elliptic integrals of Eqs. 181 to 183, the 
following expressions for the integrals are obtained: 



L 0 = A A [F(A,,<7) - F(h,q)l Ci<A (198) 



Note that ¥(X 2 ,q) = 0 if £> = A. 

L 2 = A \E(X\,q) - E(X 2 ,q)], ( 2 <A (199) 

Again, E(h,q) = 0 if ( 2 = A. 

L 4 =AI2 [2(A 2 + B 2 )E(X I ,q)-B 2 ¥(X ! ,q)] + C ! /3 [(A 2 - C, 2 )(C, 2 - B 2 )] m - 
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A/3 [2(A 2 + B 2 ) V{X 2 ,q) - B 2 F(X 2 ,q)] + C 2 /3 [(A 2 - C/)(6 2 - B 2 )] u \ Z 2 <A 



(200) 



Note that the entire second line is zero if C2 = A. 



P = 



1 



A(a 2 - A 2 ) 



n 



' A 2 -B 2 



A 2 -a] 



' J 2 A2 ~ B2 ^ 
-n A 2> — -,q 



J 



A 2 -a 



(201) 



The second elliptic integral with X 2 as its first argument is zero if C2 = A. 
[00388] Lastly, 



E = A/3 [(A 2 + B 2 ) E(X,,q) - IB 2 F(X,,q)] - {j/3 [(A 2 - C, 2 )<£, 2 - B 2 )] m - 
A/3 [(A 2 + B 2 ) E(X 2 ,q) - IB 2 F(X 2 ,q)] - C 2 /3 [(A 2 - - B 2 )) 1/2 , ( 2 <A 



(202) 



Again, the entire second line of this expression is zero if £2 = A. 

Explicit solution for interior volume V in when a s <a c + p s 

9 7 7 77 

[00389] In the case where a s < (a c + p s ) , A is still defined as it was before, A = a s 
iflc - Ps) 2 , but B 2 becomes the negative of its previous definition: 



B 2 = (a c +p s ) 2 -a s 2 . 



(203) 



7 7 

[00390] Although A is always greater than zero, we simply have that B > 0. 
The radical is defined differently as 



R V2 = j(A 2 -£ 2 )(C 2 + B 2 ) 



(204) 



and that changes the elliptic integral solutions of L 0 , L 2 , L 4 , P, and E. Defining 



Y\ = sm " 



,A p 2 + tf 
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(205a) 



y 2 = sin ' 



£i l A2 + B2 

A \B 2 + C 



(205b) 



2 J 



r = 



yjA^TB 2 



(205c) 



the following expressions are obtained for the integrals. Note that Ci is inside absolute value 
signs so that only its magnitude is used. The upper limit £?, however, may be positive, negative, 
or zero. A negative results in a change in sign for all terms in the integrals involving Q 2 - All 
the integrands are even in C so that all the integrals are odd. The elliptic integrals are also odd in 
their first argument. The expressions for the desired integrals are now 



1 



■[F(y u r) + F(y 2t r)] 



(206) 



B 2 



L 2 = 4A 2 + B 2 E(y u r)- _ — 

yj A 2 + B 

JA 2 + B 2 E(y 2 ,r)- B * 



F(mr)-\Cl 



B 2 + £ 2 



(207) 



L < = , ll D 2 K 2 * 2 - A2 ) B2F ^,r) - 2(5 4 - A*)E(y lt r)] - 

3-J A + B 



^(2A 2 - B 2 + ) 



+ 



1 



A 2 + B : 



■ [(2B 2 - A 2 )B 2 F{y 2 , r) - 2(B< - A 4 )E(y 2 , r)] - 



3 



H(2A 2 -B 2 + C 2 ) 



<B 2 + g 



(208) 
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P = 2,2 ,2,1,2 n 2 "> r ) + «. 2 ^(r ., r) + ^ 2 n(y 2 , », r) + a] F(r 2> r)] (209a) 

a 2 (a s 2 + B 2 ) y fA Y +¥ 



where 



and 



A\a] + 5 2 ) , 
« = — 7- (209b) 



E = ^A' + B 2 [B 2 F( ri ,r)-(B 2 - A>)E{y x ,r)} + ^ (£ 2 + 25 2 - ^ 2 ).f|^ 
+ i^W^,') - (2? 2 - A>)E(y 2 ,r)] + + 25 2 - A\\^L 



(210) 



which completes the solutions for the intersection of a sphere and cylinder. The explicit 
expressions for Ci and C2 are given earlier. These expressions are used in the preferred 
implementation to determine the intersection. 

[00391] The task of determining the packing fraction also may be accomplished 
according to an aspect of the invention using a modified approach. 

[00392] Recall that the radial distance of the sphere center from the z axis by 



Ps = [xs l +ys 2 ] m . (211) 

[00393] Also recall that there are four possible cases for this overlap problem. In 
Case 1, the sphere is completely outside the cylinder so V in = 0. In Case 2, the sphere is 
completely inside the cylinder so V in = (4/3)na s 3 . In Case 3, the cylinder is completely inside the 
sphere so V in = nac 2 (z t - z b ). In Case 4, the sphere and cylinder intersect each other and the 
expression for V in is more complex as will be derived below. 
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Case 1 : Sphere Completely Outside Cylinder 
[00394] The sphere is completely outside the cylinder 

V in = 0 (212) 

if one of the following conditions are met: 

If the sphere's north pole is lower than the cylinder bottom . . . 

z s + a s <zt (213a) 

Or if the sphere's south pole is above the cylinder top . . . 

z s ~a s >z t (213b) 

Or if the sphere is more than its radius away from the cylinder side . . . 

p s -a s >a c (213c) 

Or if the sphere center is above the cylinder top and is more than its radius away 
from the upper cylinder edge . . . 



2 2 2 

(Ps - o-c) + (z s - z t ) > a s and z s > z t and p s > a c (213d) 



Or if the sphere center is below the cylinder bottom and is more than its radius 
away from the lower cylinder edge . . . 



(Ps-a c f + {z b -z s f> a s 2 and z s <z b and p s >a c . (213e) 
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Case 2: Sphere Completely Inside Cylinder 
[00395] The sphere is completely inside the cylinder 

V in = (4/3)na s 3 

only if all three of the following conditions are simultaneously met: 

If the sphere's north pole is no higher than the cylinder top . . . 

z s + a s < z t 

And if the sphere's south pole is no lower than the cylinder bottom 

z s -a s > Zb 

And if the sphere's equator nowhere exceeds the cylinder radius . . 

p s + a s <a c . 

Case 3: Cylinder Completely Inside Sphere 
[00396] The cylinder is completely inside the sphere 

V in = na c 2 (z t - z b ) 

if both the following two conditions are met: 

If the upper edge of the cylinder is inside the sphere . . . 

(Ps + <*c) 2 + (zt-Zs) 2 <ci s 2 
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And if the lower edge of the cylinder is inside the sphere . 



(p s + a c ) 2 + (z s -z b f<a s 2 . 



(217b) 



Case 4: Cylinder and Sphere Intersect 

[00397] Again, there are several subsets of Case 4. To efficiently organize them, we 
define two circles: the sphere circle and the cylinder circle. The sphere circle is the perimeter of 
the sphere's cross section at a given z altitude. Its radius is 



and its radial position p s is given by Eq. 21 1. Likewise the cylinder circle is the perimeter of the 
cylinder's cross section which has constant radius a c at any z altitude between z b and z t and is 
centered on the z axis. 

[00398] To describe the intersection of the sphere and cylinder, we define four subcases 
for the relationship of the sphere circle and the cylinder circle at a particular altitude z: 



[00399] To find V in , our task is reduced to integrating in z the overlap areas of the two 
circles A in from the altitude where one subcase holds to the altitude where another subcase holds. 



r s = [a 2 - (z - z s f] m , z s -a s <z<z s + a s 



(218) 



Subcase 4.1 : The sphere circle is completely enclosed in the cylinder circle. 
Subcase 4.2: The two circles do not overlap or one or both are nonexistent. 
Subcase 4.3: The cylinder circle is completely enclosed in the sphere circle. 
Subcase 4.4: The two circles intersect. 



For Subcase AA,A i} 



in 




(219a) 



For Subcase 4.2, A it 



in 



0. 



(219b) 
(219c) 
(219d) 



For Subcase 4.3, A it 
For Subcase 4.4, A u 



in 



in 



na 2 . 

2 2 

a c {a - cosa sina) + r s {fi - cos/3 sin/?) 



where r s is given by Eq. 218 and 
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cos a = 



Ps +<*] ~rl 



(220a) 



sin a = 



Vfc 2 -(a e -p t ) 2 ][(a e + p t ) 2 -rn 
2p s a c 



(220b) 



cos P = 



Ps + r s -<*t 
2p s r s 



(221a) 



sin P - 



4\r] -(a c - Ps y][(g c + Ps Y-rn 
2Ps r * 



(221b) 



[00400] It is now necessary to determine a series of altitudes at which the cylinder circle 
and the sphere circle change from one subcase to another in Case 4. 

[00401] A sphere and a cylinder side can intersect each other at either 0, 2, or 4 altitudes. 
Pairs of altitudes may be degenerate (two solutions having the same altitude) if the sphere just 
nicks the cylinder wall. The solutions for these altitudes come from solving the following 
equation: 



[r s (z)] 2 = [a c ±p s f 



(222a) 



or 



a s 2 -(z-z s f = [a c ±p s ] 2 



(222b) 



which, using Eq. 218, has the four solutions 




(223) 
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The + sign in the first ± of Eq. 223 denotes the upper or northern hemisphere intersection 
altitude; the - sign denotes the lower or southern hemisphere intersection altitude. The + sign in 
the second ± of Eq. 223 denotes intersection with the far cylinder wall, on the opposite side of 
the z axis from the sphere's side; the - sign denotes intersection with the near cylinder wall, on 
the same side of the z axis as the sphere. 

[00402] Armed with these four solutions, we will now consider each of the four 
subcases of Case 4, but we will consider separately whether the sphere circle lies in the southern 
hemisphere of the sphere or whether it is the equatorial plane or northern hemisphere. We 
consider the four subcases separately for the two hemispheres making eight situations in all. We 
consider the hemispheres separately because if the sphere circle is in the southern hemisphere, it 
is increasing in size as z increases. Hence the preferred implementation can determine the 
altitude at which subcases change. If the sphere circle is its equatorial plane or is its northern 
hemisphere, it is decreasing in size as z decreases, and the altitude at which subcases will change 
again can be determined. 

[00403] We now outline the overall logic used in the preferred implementation to find 
the altitudes at which the subcases for Case 4 change. Let z/ ovv be the lower altitude of the 
subcase region currently being examined and Zhigh be its upper altitude. 

[00404] The absolute boundary below or above, with never any overlap between sphere 
and cylinder, is given by 

z 0 = MAX[z b , z s - a s ] (224) 

for the lower boundary and 

z 5 = MIN[z„ z S9 + a s ] (225) 

for the upper boundary. The altitudes z\ through Z4 correspond to the different intersection 
altitudes (if they exist) of the sphere's polar plane that is perpendicular to the cylinder axis. A 
sphere's polar plane is the plane which intersects the sphere in a great circle perpendicular to the 
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equatorial plane. It can be oriented at any azimuthal angle. The polar plane of interest to us is 
the one passing through the cylinder axis. 

[00405] In ascending order (again, if they exist) the intersection altitudes are: 



Southern hemisphere intersection with cylinder wall closest to sphere center: 



*i = W«? ~ ( a c ~ Ps) 2 (226) 
Southern hemisphere intersection with cylinder wall furthest from sphere center: 
z 2 = z s -^a 2 s -{a c + p s f (227) 
Northern hemisphere intersection with cylinder wall furthest from sphere center: 
z 3 = *, + V«, 2 -{a c + p s f (228) 
Northern hemisphere intersection with cylinder wall closest to sphere center: 
z A = z s + 4a] -(a c -p s f . (229) 

In all four of these solutions, if the argument of the square root is negative, the square root is 
replaced by zero. If the sphere intersects both sides of the cylinder, there are four real solutions; 
if only one side, there are two real solutions; if not at all, there are no real solutions of Eqs. 226 
through 229. The figure below shows the two real solutions for a sphere that intersects only one 
side of the cylinder. 

[00406] The method used in the preferred implementation for finding the volume V in of 
a sphere that lies inside a cylinder will now be discussed. The method is divided into two DO 
WHILE loops, one for the southern hemisphere and one for the northern hemisphere (which 
includes the sphere's equatorial plane). This is convenient because for the southern hemisphere 
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the sphere circle is always increasing with increasing z while for the northern hemisphere it is 
always decreasing. The method as specifically implemented here has a number of comments so 
that the reader can follow the reasoning guiding it. 



Vin = 0 ! Vin is the volume of this sphere lying inside the 

! cylinder. Initialize it to zero, 
zo = MAX[Zb> z s - as] ! zo and z 5 are the absolute lower and upper 

Z5 = MIN[z t , z s + as]! limits of possible sphere-cylinder overlap. 



IF (z s > Zb) THEN ! The southern hemisphere can contribute to Vj n if and 

! only if this statement is true. Otherwise it lies 

! beneath the cylinder bottom and cannot contribute. 

! (It might not contribute even when the above 

! statement is true but that is checked below.) 



ziow = zo ! Initialize zi ow to smallest possible value. 

r s (ziow) = -(ziow-Zs) ] ! Find initial sphere circle 

! radius. It is 0 if square 
! root argument < 0. 

! Find the initial southern hemisphere subcase. 
IF (p s + r s (zi ow ) < ac) THEN 

subcase = 1 ! Sphere circle inside cylinder circle. 
ELSE IF (p s - r s (z low ) > ac) THEN 

subcase = 2 ! Sphere circle outside cylinder circle. 
ELSE IF (r s (ziow) - Ps > ^) THEN 

subcase = 3 ! Cylinder circle inside sphere circle. 

ELSE 

subcase = 4 ! Sphere circle intersects cylinder circle. 
END IF 
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! We exit the DO WHILE loop when the subcase is no longer a positive 
! integer. Remember that since we are in the southern hemisphere, the 
! sphere circles are expanding with increasing z. 

DO WHILE (subcase > 0) ! Southern hemisphere loop. 

IF (subcase = 1) THEN 

! The sphere circle lies inside the cylinder circle. 
! It expands with increasing z until one of three 
! things happen: 

! (1) It hits the top of the cylinder. 

! (2) It intersects the cylinder circle. 

! (3) It reaches its equatorial plane. 

! Numbers (2) and (3) can be combined if we agree to 

! define z\ with the square root replaced by zero if 

! its argument is negative. 

w = as 2 - (ac - p s ) 2 
IF (w > 0) THEN 



Zi =z s 
END IF 

z high = MIN[zt, zi] 

Vj n — Vj n 

+ Vi(zi 0WJ Zhigh) ! We'll derive an explicit 
! expression for this Vi later. 

! Now we must find the next subcase. 

IF (z high < MIN[Zs,zJ) THEN ! There is another 

! subcase to be 
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IF (p s > 0) THEN 
subcase = 4 



ELSE 



subcase = 3 



END IF 



ELSE 



subcase = 0 



END IF 



! addressed in the 
! southern hemisphere. 
Sphere is not concentric with 
the cylinder so it will 
expand until it intersects. 
Here, the sphere is concen- 
tric so it will expand 
until it perfectly circum- 
scribes the cylinder 
circle. 

! You reached the end of southern 
! hemisphere contribution. End the 
! S. hemisphere's DO WHILE loop. 



ELSE IF (subcase = 2) THEN 

The sphere circle lies outside the cylinder circle 
until it expands to the point where it intersects 
the cylinder. This intersection must happen or we 
wouldn't be here in Case 4. The intersection 
happens at z\. There is no contribution to Vi n from 
this region. We only need to find Zhi g h = z\. 



w = as 2 - (ac - ps) 2 
IF (w > 0) THEN 



,1/2 



Zhigh = z s - w 

subcase = 4 ! Not done yet with S. hemisphere. 



ELSE 



Zhigh Z s 

subcase = 0 ! Reached equatorial plane. Now 
END IF Iwe're done with S. hemisphere. 
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ELSE IF (subcase = 3) THEN 

! The cylinder circle lies inside the sphere circle. 
! Since the sphere circle is expanding, it will 
! remain inside until it reaches the first of either 
! the top of the cylinder or its equatorial plane. 

z high = MIN[zt, z s ] 

Vin = Vj n + V2(zi OW5 Zhigh) ! We'll derive an explicit 

! expression for this V2 later. 

subcase = 0 

ELSE EF (subcase = 4) THEN 

! The only possibility left is that sphere and 

! cylinder circles initially intersect. This region 

! extends upward until one of three things happen: 

! (1) We reach the top of the cylinder. 

! (2) We reach the sphere's equatorial plane. 

! (3) We reach a point at which the sphere circle 

! completely circumscribes the cylinder circle. 

! We cannot reach the top of the sphere in this case 

! because we are only dealing with the southern 

! hemisphere. (Remember Z2 = z s if the argument of 

! the square root is negative.) 

w = as 2 - (ac + p s ) 2 
IF (w > 0) THEN 

z 2 = z s - w 1/2 

subcase = 3 

ELSE 
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z 2 = z s 
subcase = 0 
END IF 

z high = MIN[zt, z s , z 2 ] 

Vj n = Vi n + V3(ziow,Zhigh) ! We'll derive a numerical 

! integration for this V3 later. 

END IF 



If you have reached either the cylinder top or the 
sphere's equatorial plane, you are done with the southern 
hemisphere's contribution to Vi n . Go to the northern 
hemisphere next. (We use > signs below in case of round- 
off error. With infinite precision, Zhigh could never be 
greater than z t or z s .) 



IF (zhigh > z t OR Zhi g h > z s ) THEN ! Redundant check saying 
subcase = 0 ! we're done with the southern hemisphere. 

ELSE 

ziow = Zhigh ! Cycle back through the DO WHILE loop 
END IF ! for the next region in the southern 

! hemisphere contributing to Vj n . 
END DO ! End of southern hemisphere DO WHILE loop. 



END IF ! End of southern hemisphere overall IF statement. 



IF (z s < z t ) THEN ! The northern hemisphere can contribute to V in if and 

! only if this statement is true. Otherwise it lies 

! above the cylinder top and cannot contribute. (It 

! might not contribute even when the above statement 

! is true.) 
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ziow = MAX[z b , z s ] 

r s (ziow) = [as 2 -(zi ow -z s ) 2 ] 1/2 ! Find initial sphere circle radius in 

! N. hemisphere. It is 0 if square root 

! argument <0. 



! Find the initial northern hemisphere subcase. 
IF (p s + r s (z ]ow ) < ac) THEN 

subcase = 1 ! Sphere circle inside cylinder circle. 
ELSE IF (p s - r s (zi ow ) > ac) THEN 

subcase = 0 ! Sphere circle outside cylinder circle. 

! Although this is actually subcase 2, since 
! N. hemisphere sphere circle can only contract 
! it will never intersect the cylinder 
! circle if it does not initially intersect it. 
! Therefore, we set the subcase to zero so the 
! N. hemisphere's DO WHILE loop will end. 
ELSE IF (p s + ac < r s (z low )) THEN 

subcase = 3 ! Cylinder circle inside sphere circle. 

ELSE 

subcase = 4 ! Sphere circle intersects cylinder circle. 
END IF 



DO WHILE (subcase > 0) ! Begin northern hemisphere loop. 

IF (subcase = 1) THEN 

! The sphere circle lies inside the cylinder circle. 
! Since the sphere circle shrinks as z increases, it 
! cannot ever intersect the cylinder circle. It can 
! only end at its north pole or at the top of the 
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! cylinder. That ending position is z 5 . 

Zhigh = z 5 

Vin = Vin + Vi(zi 0W3 Zhigh) ! We'll derive an explicit 

! expression for Vi later. 

subcase = 0 

ELSE IF (subcase = 2) THEN 

! The sphere circle is inside the cylinder circle. 
! In the northern hemisphere this means they will 
! never intersect more so there will be no more 
! contribution to Vi„. 

subcase = 0 

ELSE IF (subcase - 3) THEN 

! The cylinder circle lies inside the sphere circle. 
! The sphere circle will contract until we hit top 
! of the cylinder or until we reach z 3 which is the 

! point at which the two circles will intersect. 

2 2 

w = as - (ac + p s ) ! w is non-negative or we 

! couldn't be here. 

z 3 = z s + w 1/2 
Zhigh = MIN[z t , z 3 ] 

Vin = Vj n + V2(ziow>z n igh) ! We'll derive an explicit 

! expression for V2 later. 

IF(z hign <z t ) THEN 

IF (p s > 0) THEN ! Sphere is not concentric 
subcase = 4 ! with cylinder so it will 



117 



contract until it inter- 



ELSE 



subcase = 1 



sects the cylinder circle. 
The sphere is concentric so 
it will contract until it 



ELSE 



END IF 



! suddenly lies completely 

! inside the cylinder circle. 
! If we're here, we reached 



subcase = 0 



top of the cylinder; we're 
done with N. hemisphere. 



END IF 



ELSE IF (subcase = 4) THEN 

! The only possibility left is that sphere and 
! cylinder circles initially intersect. This region 
! extends upward until one of two things happen: 
! (1) We reach the top of the cylinder before 
! reaching the north pole. 
! (2) The intersection stops with the sphere circle 
! now either completely inside or completely 
! outside the cylinder circle. This happens at 
! z 4 . 

w = as - (ac - p s ) ! w is never negative or we 



z 4 = z s + w 

Zhigh = MIN[z t , z 4 ] 

Vin = V in + V 3 (zi ow ,Zhi g h) ! We'll derive a numerical 



IF (p s < ac) THEN 

subcase = 1 ! Next, sphere circle will lie 



wouldn't be here. 



solution for V 3 later. 



! inside cylinder circle. 
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ELSE 

subcase = 0 ! Next, sphere circle will lie 
! outside cylinder circle. 

END IF 
END IF 

! (Note the second logical expression in the EF statement below is zs, 
! notz s .) 

IF (zhigh > Zt OR Zhigh > Z5) THEN ! Redundant check saying 
subcase = 0 ! we're done with northern hemisphere. 

ELSE 

ziow = Zhigh ! Cycle back through the DO WHILE loop 
END IF ! for the next region in the northern 
! hemisphere's contribution to Vi n . 



END DO ! End of northern hemisphere DO WHILE loop. 
END EF ! End of northern hemisphere overall EF statement. 

[00407] Expressions for the three volume expressions F/(z/ ow , Zhigh), i = 1, 2, or 3 will 
now be addressed. Remember in the following that a c is the cylinder radius, a s is the sphere 

radius, (x s , y s , z s ) is the position of the sphere's center, and p s = y]x* + y] is the radial distance 

of the sphere center from the cylinder axis. 



DERIVATION OF V } 

[00408] The volume V\ is an integral of full sphere circles from z/ ow to z^gh because the 
sphere circles lie entirely within the cylinder circle. 



z high 



V x = n J dz[r s (z)f = n J dz[a] - (z - z s ) 2 ] (230) 



z low 
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or 




4* + z«g* (3z s - z low ) - 3z s (Zj - z to J] 



(231) 



where z/ 0H , and Zhi g h depend on where we are in the above method. Their values are given in the 
implementation shown above. 

DERIVATION OF V 2 

[00409] The volume Vi is the integral of the cylinder circle from z/ ovv to z^ because the 
cylinder circle lies entirely within the sphere circle. It is given simply as 



DERIVATION OF V 3 

[00410] The volume V3 is the integral of intersecting sphere and cylinder circles from 
zi ow to Zhigh The limits are always given in the method as described above. The intersect area is 
given by Eqs. 219, 220, and 221 which we must integrate in z. 



V 2 = 7ta 2 c (z high - z low ) . 



(232) 



F 3 = J dz{a\ [a - cos a sin a] + [a 2 s - (z - z s ) 2 ][/3 - cos ft sin /?]} 



(233) 



'law 



where 



cos a = 



2 Ps°c 



(234) 



sin a = 



-{a c -p s y-(.z-z s y][{a c + p s y-al +(z-z s y] 

2 Pfic 



(235) 
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a = arccos 



A 2 +" 2 c -at +(z-z s Y 



(236) 



and 



cos P = 



A 2 - a 2 +a 2 5 



-{z-z s Y 



Ip^a] -{z-z 5 f 



(237) 



sin /3 



V[g 5 2 - (g c - a) 2 - (z - z s f][{a c + A) 2 + « 2 - (* - 
2 A A 2 -(^-^) 2 



(238) 



/? = arccos 



A 2 -« c 2 +a 2 



- (* - z 5 ) 2 



2p^a\ -(z-z s f 



(239) 



[00411] The integral of Eq. 233 could be done analytically in terms of complete and 
incomplete elliptic integrals of the first, second, and third kinds. However, this solution is so 
cumbersome that it is more efficient to perform it numerically, especially because the integrand 
is quite smooth and well-behaved so that a Romberg integration is very efficient and contains an 
accurate error assessment. 

Calculating the Total Volume Fraction for a Pack 

[00412] The preferred implementation also can determine the volume fraction § of the 
pack. Only the innermost cylinder contains particles that look like a full three-dimensional pack. 
The outer cylinders do not contain the particles the inner cylinders contain. However, it would 
be unfortunate to throw away all the pack information the outer cylinders contain. They do 
perturb the pack of the innermost cylinder correctly and their correct geometry for the particles 
they contain should somehow be used in estimating the volume fraction of the pack. 

[00413] The preferred implementation determines the volume fraction ^ in the 
following way. The y'-th cylinder can contain all particle modes whose cylinder radii are less 

than or equal to its cylinder radius cc c (/). As long as the water level option of pack building is 

121 



operating, the y-th mode particles are distributed approximately as they would be in a truly 
three-dimensional pack. The preferred implementation therefore calculates the volume fraction 
$j of Mode j particles within the y-th cylinder and calls it the y-th partial volume fraction. It does 

likewise for all the other modes, e.g., J of them. Then the total volume fraction is simply the 
sum of all the partial volume fractions: 



If desired, the user can specify that only the volume fraction of the innermost cylinder, where the 
pack is three dimensional, will be calculated. 

[00414] To calculate the partial volume fraction <f> } in the y-th cylinder, the preferred 

implementation sums the interior volume of each Mode j sphere inside the y-th cylinder. Let nj 
be the number of Mode j spheres in the y-th cylinder. Let Vj(k) be the volume inside the cylinder 
of the k-th sphere of Modey. Then 



[00415] Dividing this total interior partial volume by the volume of the y-th cylinder 
Vcyiij) gives the partial volume fraction <j) } for the pack: 



j 



(240) 



(241) 



*=1 



h = V lot <J)/V cyl (j) 



(242) 



where 



V C yi(j) = Tt[a c (j)] 2 [z, op (J) - z bot (j)]. 



(243) 



[00416] The user has the option of making calculating volume fractions in cylinders that 
are a little smaller than the actual cylinders used in making the pack so that possible wall effects 
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can be eliminated from the calculation. Remember, however, that the cylinder walls are barriers 
to the particle centers, not their leading edges so that disturbances to the pack from wall effects 
normally will be small. 

Finding the Radial Distribution Functions 

[00417] The radial distribution functions, denoted by g/,(r)dr, are probability distribution 
functions denoting the mean number of Mode or Submode j particles whose centers are between 
r and r + dr from the centers of Mode or Submode i particles throughout the pack. They can be 
important in establishing many properties of random packs. Because the packs here are assumed 
random, only the scalar radial distance r is used rather than vector r in this implementation. 

[00418] Reduced dimension packs do not extend equally in all directions. They are a 
series of concentric cylinders. Therefore, for a given Mode or Submode i particle, the preferred 
implementation calculates the volume of the spherical shell segment centered on Particle / and 
lying between r and r + dr, that lies inside the y-th cylinder. Then it counts all Mode j particles 
lying in that shell segment. Dividing the number of Mode j particles by the volume of the shell 
segment gives an approximation to gij(r)dr at that position r. 

[00419] Finding the volume of this spherical shell segment is not a trivial procedure. 
Fortunately, a means to calculate the total sphere volume lying inside any cylinder is described 
above. Therefore, the preferred implementation can calculate the total volume of an imaginary 
sphere of radius r + dr, centered at Particle /, that lies inside the y-th cylinder, and subtract from it 
the total volume of the concentric sphere of radius r to obtain the volume of the spherical shell 
segment between r and r + dr lying in the y-th cylinder. 

[00420] The user specifies the maximum radial distance r max (ij) for which gy is desired 
and the number of bins nty into which it divided r max (ij). Then the bin widths have finite width 

Ar, y = r max (ij)lmij (244) 

rather than infinitesimal width dr. Then the gy are determined as follows. 

[00421] First, note that the g,y(r) are zero for r < a t + aj. Then, consider a particular 
Mode /' particle with index number //,. Denote the number of Mode j particles in the first bin 
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from Particle ih particle by NJJh\ «/ + aj). The first bin lies between (a,- + aj) and (a, + aj) + Ar ; y. 
In general, the number of Mode j particles lying in the A>th bin from Particle ih is denoted by 
Nj(ih\ fk) where n = ai + aj + (fc-l)Ar,y, £=1,2,..., 

[00422] Now consider the volume of the &-th bin from Particle ih that lies inside the y'-th 
cylinder. This can be obtained from the volumes V in of previous sections where Vi n = V in (x s , y s , 
z s , a s \ a c , Zb oh z top ). The argument list contains the position of the sphere center (x S9 y s , z 5 ), the 
sphere radius a s , the cylinder radius a C5 and the cylinder's bottom and top altitudes {z^ oU z top ). 
The volume for the k-th bin about Particle ih is then given by 

Vin(ihj\ k) = V in {xih, yjh> Zih, n+h Gcj, Zbotj, Ztopj) - V in {xih, ym, Zfa Tk\ a c j, ztotj, z top j) (245) 

and the number density of Mode j particles in this volume is 

o(ihj\ k) = Nj(i h \ r k )IV in {i h J\ k) . (246) 

[00423] Finally, to obtain the g,y(r), the above expression can be averaged over all Mode 
/ particles: 

g {j {r)~-Y j a{i h J\k) (247) 
"i in 

if r > a.i + a.j and is zero otherwise. The term k is the integer lying in the range 

(r - fly - qD/Artj < k < 1 + (r - a t - a^/Anj . (248) 

If r is so large or so small for a particular particle i h that the spherical shell containing it misses 
the y'-th cylinder altogether, then that term in the sum of Eq. 247 is omitted altogether. 

[00424] The normalized radial distribution function, denoted G,y(r), is given by 

Gij(r)=g i j{r)/(4Kr 2 ) (249) 
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or in our finite bin width treatment, it is given by 



-x(r l k + x -rl) 



(250) 



It is the mean number density nj of Type j particles at a distance r from a Type / particle. Thus 
normalized, the function Gy(r) approaches the constant nj as r grows large rather than increase 
quadratically as gy{r) does, nj represents the number of j particles per unit volume, e.g., 
concentration, output from the preferred implementation. 

Finding the Coordination Numbers 

[00425] The coordination number Qy is the mean number of Mode j particles touching a 
Mode i particle. Again, because a reduced dimension pack is ordered in concentric cylinders 
rather than being isotropic in all directions, the preferred implementation finds how much of the 
surface of the z-th particle is available for Mode j particles to touch. That is equivalent to asking 
how much of the imaginary spherical surface centered on the Mode / particle and having radius 
a,- + Oj lies inside the y'-th cylinder. We will denote that interior portion of the surface by Sy. 
Then if a search shows there are hy Mode j particles touching a particular Mode i particle 
identified by index number 4, then our best estimate of the total number of touching Mode j 
spheres, if the entire surface were accessible, would be 

W = h ff 4n(ai + aj) 2 /Sij . (251) 

[00426] If the entire surface is accessible, then Sy = 4rc(<2, + ajf and fjy = hy. Because of 
round-off error, a certain center-to-center error tolerance 8 is required to determine whether two 
spheres are indeed touching. The code considers them touching if their centers are within the 
distance (a ( - + aj)(l + e) from each other. Because the preferred implementation is in double 

* • 12 

precision, the default for e is 10" , but the user can specify any other desired 8. Normally it 
should be a few orders of magnitude greater than machine precision. 
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[00427] The spherical surface Sy lying inside a cylinder could be obtained from our 
earlier work on obtaining the spherical volume inside a cylinder by taking a derivative of the 
volume expressions with respect to the sphere radius. However, those expressions are 
complicated and have substantial bookkeeping associated with them. Because in this 
implementation there is an error tolerance, we may take a numerical derivative as follows: 

Sy = G l [Vin(Xih 9 yih, Zfa fa + aj){\ + e); a c (j\ z botJ , z topJ ) - 

Vinfah* yih Zih, (fli + Ay); a C {j\ Zbotp Ztopj)] - (252) 

The £ used in this expression should be kept very small even if the user chooses to be rather 
liberal in assigning a rather large error tolerance. One might, for example, assign this 8 to be 
either the user-defined (or default) error tolerance or 10~ 6 , whichever is smaller. 



Finding Various Pack Statistics 

[00428] There are a number of vital statistics of the pack that will generally be of 
interest. If the user so specifies, then the following data is written to a user-named file. 

[00429] 1 . A recapitulation of the particle histogram. 

[00430] 2. A recapitulation of the values for control flags and parameters 
used in building the pack. 

[00431] 3. The total number of particles placed in the pack and the number 
placed from each mode. 

[00432] 4. The minimum and maximum radial and axial positions for the 
centers of all particles in the pack and the same for each mode. 

[00433] 5. Radial and axial particle density functions for each mode 
specified between the minimum and maximum positions specified above. 

[00434] 6. The number of particles placed after tunneling (if tunneling is 

permitted). 

[00435] 7. The number of particles hitting one of the catch nets. 
[00436] 8. The number of particles hitting one of the water levels. 
[00437] 9. The number of preplaced particles, if any. 
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[00438] It is important to understanding in implementing the various aspects of the 
invention and in applying the related principles that the mathematical equations and expression 
presented herein may be assumed to have an infinite degree of precision and all real 
implementations of those equations result in answers with finite precision. Many programmers 
assume that the use of double precision, 64-bit precision, eliminates the need for close attention 
to this fact. Such an assumption with this problem typically will be incorrect. For example, in 
testing the implementation of several equations, cases were found where use of quadruple 
precision, 128-bit precision, was insufficient to obtain correct solutions. These cases were 
avoided if in the solution process checks were made that allowed for tolerance bands. Every 
check for equality must include equality within a tolerance. The following tolerances are used 
within the preferred implementation. 

[00439] 1 . Comparison with zero. Numbers in the range -1 .OE-20 < number 
< 1. OE-20 are treated as zero. 

[00440] 2. General comparison for equality. Any two non-zero numbers are 
within 1.0E-10 times the smaller of the two numbers are treated as equal. 

[00441] 3. General angle comparison for equality. Any two angles within 
1.0E-6 radians are treated as equal. 

[00442] 4. Distance comparison. Any two distances within 1.0E-6 times the 
smallest particle radius are treated as equal. 

[00443] 5. Distance squared comparison. Any two numbers representing 
distance squared within 1.0E-6 times the square of the smallest particle radius are treated as 
equal. 

[00444] Additional advantages and modifications will readily occur to those skilled in 
the art. Therefore, the invention in its broader aspects is not limited to the specific details, 
representative devices and methods, and illustrative examples shown and described. 
Accordingly, departures may be made from such details without departing from the spirit or 
scope of the general inventive concept as defined by the appended claims and their equivalents. 
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